Introduction

Interactions of brain oscillations in different frequency bands, also known as cross frequency couplings (CFC), have been proposed to serve as a mechanism for neural coordination across different spatiotemporal scales (Canolty and Knight, 2010). The most common form of CFC is phase-amplitude coupling, whereby the phase of the lower frequency rhythms modulates the amplitude of spatially localized high frequency oscillations. This type of modulation allows for efficient control of communication and information transfer between brain regions. In humans, phase-amplitude CFC has been observed in multiple brain areas based on recordings obtained at different levels, i.e., from intracellular to surface electroencephalogram (EEG) measurements and under different experimental conditions (Mormann et al., 2005; Canolty et al., 2006; Cohen et al., 2008; De Lange et al., 2008; Penny et al., 2008; Tort et al., 2009; Axmacher et al., 2010; Seeber et al., 2014; Gwon and Ahn, 2021) or pathological conditions (López-Azcárate et al., 2010; Özkurt and Schnitzler, 2011; De Hemptinne et al., 2013; Edakawa et al., 2016; Combrisson et al., 2017; Wang et al., 2017).

In non-invasive brain computer interfaces (BCI), the concept of CFC has not yet been fully explored. To date, only a few studies have investigated the prospect of CFC in BCI systems (Gysels and Celka, 2004; Wei et al., 2007; Dimitriadis and Marimpis, 2018; Georgiadis et al., 2019; Feng et al., 2020; Gwon and Ahn, 2021). Gysels and Celka (2004) and Wei et al. (2007) used CFC as a measure of functional connectivity between different brain regions and demonstrated its potentiality as a feature for mental task classification. In Dimitriadis and Marimpis (2018), CFC was estimated within EEG sensors from single trials during experimental paradigms that generated visually evoked potentials. By employing CFC as a feature for decoding flashing images and stimulus presentations, high classification accuracy and information transfer rates were achieved. Similarly, in Georgiadis et al. (2019), Feng et al. (2020), and Gwon and Ahn (2021) CFC was introduced as the main building block of motor imagery based decoding schemes.

When it comes to CFC estimation, there is no golden rule that applies in selecting an appropriate method. Some of the most commonly used CFC estimation techniques are the Phase-locking Value (PLV) (Lachaux et al., 1999; Cohen, 2008), the Mean Vector Length (MVL) (Canolty et al., 2006), the Kullback–Liebler based Modulation Index (MI) (Tort et al., 2010), and the Generalized Linear Model (GLM) (Penny et al., 2008). Each method was developed based on different principles and exhibits unique advantages but also limitations. The PLV compares the phase of the instantaneous high frequency amplitude with the low frequency phase. The MVL computes the magnitude of the average complex signal defined by the instantaneous high frequency amplitude and the instantaneous low frequency phase. The MI measures the deviation of the distribution of the high frequency amplitude, with respect to the phase of the low frequency, from the uniform distribution and the GLM is simply a linear regression model fitted on the instantaneous high frequency amplitude using the cosine and the sine of the instantaneous low frequency phase as regressors.

Herein, we provide an alternative CFC estimation method based on Linear Parameter Varying Autoregressive (LPV-AR) models (Tóth, 2010). CFC is a non-linear phenomenon and therefore requires non-linear methods for its estimation. In contrast to the simple linear AR models, LPV-AR models are able to capture non-linear dynamics and interactions. The unique characteristic of these models is that their coefficients can be modulated by an external variable/signal and, more importantly, they can be identified using the Least Squares approach since the LPV-AR is linear in the parameters. The main idea is to use the instantaneous low frequency phase as the external modulating signal of an AR model fitted on the instantaneous high frequency amplitude. AR models have been widely employed for power spectrum estimation. Specifically, the AR coefficients can be used to describe the power spectrum of a signal. By allowing an external variable to modulate these coefficients, we are modeling the time-varying spectral changes of the signal under consideration induced by the external variable. There has been one previous study on LPV-AR based CFC (la Tour et al., 2017), however, their approach in estimating CFC is completely different than the one that we describe here. The model of La Tour et al. is applied directly on the raw signals, in contrast to the instantaneous amplitude and phase as usually done in CFC estimation. However, filtering and preprocessing is required in order to isolate the low-frequency driver and remove its effects from the raw signal. The model also assumes time-varying model innovations which requires the use of non-linear optimization techniques for model estimation. On the other hand, our proposed approach provides closed-form solutions to a linear Least Squares problem.

Our proposed LPV-AR based CFC methodology was validated using both synthetic and real data. The real data comprised of EEG recordings obtained during hand/arm movement attempts of spinal cord injured (SCI) participants. Our main goal is to encourage the use of CFC as a decoding feature in future BCI applications but also validate the capabilities of our proposed LPV-AR based CFC measure using real EEG recordings.

Materials and methods Linear parameter varying autoregressive model based cross frequency coupling

In this section we describe the LPV-AR modeling technique and its application on CFC estimation.

Linear parameter varying autoregressive model

In an autoregressive (AR) model, the signal of interest is expressed as a linear combination of its past values (Ljung, 1998),

y⁢(n)=∑k=1pak⁢y⁢(n-k)+ε⁢(n)(1)

where n represents the discrete time, y ∈ RN×1 is the observed signal, N is the total number of samples, p is the model order (i.e., the number of past time lags that are taken into account), ak are the AR coefficients for each order k and ε ∈ RN×1 is a zero-mean white noise vector. In LPV-AR systems, the coefficients ak are modulated by a time-varying external signal s ∈ RN×1 known also as scheduling variable. Eq. 1 is then expressed as (Tóth, 2010),

y⁢(n)=∑k=1pak⁢(sn)⁢y⁢(n-k)+ε⁢(n)(2)

where for notational simplicity sn = s(n). Note here that ak has a static dependence on s (i.e., it only depends on the instantaneous values of s at each time point n). Identification of the LPV-AR model of Eq. 2 can be addressed by recasting the problem as linear regression. Specifically, the AR model coefficients are expressed as a linear combination of a set of known fixed basis functions,

ak(sn)=θkTΨk(sn)=θk0+∑i=1qθkiψki(sn),k=1:p(3)

where q is the number of basis functions. Herein, we selected a polynomial basis due to its simplicity. Compared to other basis functions (e.g., radial basis functions) that require tuning of multiple hyperparameters, the polynomial basis requires only the selection of a polynomial order q, e.g., ΨkT⁢(sn)=[1,sn,…,snq]. Furthermore, polynomial functions are widely used in biosignal processing and physiological systems (Mitsis and Marmarelis, 2002; Marmarelis, 2004; Zhang et al., 2010). Based on the polynomial basis, Eq. 2 can be expressed as,

y⁢(n)=wT⁢φ⁢(n)+ε⁢(n)(4)

where,

w=[θ10,θ11,…⁢θ1q,…,θp0,…⁢θpq]T∈RD×1(5)

with D being the total number of coefficients i.e., D = p(1 + q) and φ(n) is the extended regressor,

φ(n)=[y(n-1),sny(n-1),…snqy(n-1),…,

y(n-p),sny(n-p)…,snqy(n-p)]T(6)

For each order p one could define a different q value, however, for simplicity we kept q fixed. If N is the total number of samples (i.e., n = 1…N), the data can be reorganized in the following form,

Φ=[φ⁢(1),…⁢φ⁢(N)]T∈RN×D(7)

In matrix form, Eq. 4 can be written as,

y=Φ⁢w+ε(8)

where y,ε ∈ RN×1. The unknown coefficients w can be estimated using the Least Squares (LS) solution,

wLS=(ΦT⁢Φ)-1⁢ΦT⁢y(9)

In addition, if we apply ridge regression, the regularized LS (RLS) solution is given as,

wRLS=(ΦT⁢Φ+λ⁢I)-1⁢ΦT⁢y(10)

where λ is the so-called regularization parameter.

Linear parameter varying autoregressive based cross frequency coupling

To quantify the CFC between the phase of a low-frequency rhythm (i.e., phase frequency band L with a frequency range of [LminLmax]) and the amplitude of a high-frequency rhythm (i.e., amplitude frequency band H with a frequency range of [HminHmax]), we first apply band-pass filtering in conjunction with the Hilbert transform (Papoulis and Pillai, 2002). Specifically, the raw signal x ∈ RN×1 is band-pass filtered in the low and high frequency bands of interest (i.e., L and H) using a two-way least-squares finite impulse response filter [Matlab function eegfilt.m from EEGLAB toolbox (Delorme and Makeig, 2004)], to preserve phase information. The filter order depends on the low-frequency cut-off value and is defined as the rounded value of three times the ratio of the sampling rate to the low frequency cut-off or else order of three cycles of the lower cut-off frequency. The obtained signals are xL ∈ RN×1 and xH ∈ RN×1. Next, we employ the Hilbert transform to create the analytic signals of xL and xH and extract the instantaneous phase PL ∈ RN×1 and amplitude AH ∈ RN×1 of the band-pass filtered signals in the frequency bands L and H, respectively,

PL=∡⁢H⁢i⁢l⁢b⁢e⁢r⁢t⁢ a⁢n⁢d AH=||H⁢i⁢l⁢b⁢e⁢r⁢t⁢||(11)

Note that the Hilbert transform expresses real-valued signals as complex functions with time-varying amplitude and phase (also known as analytic signals). The analytic signals are useful in extracting instantaneous attributes of different time-series. The resulting amplitude and phase envelopes of Eq. 11 are used to assess the phase-amplitude interactions between the low and high frequency bands of interest. To quantify the CFC between the phase of the L band and the amplitude of the H band, we use the LPV-AR methodology described in section “Linear parameter varying autoregressive model.” The scheduling variable s is assumed to be the cosine and the sine of the instantaneous phase in the L band, whereas the signal of interest y is the instantaneous amplitude in the H band normalized by its norm (to remove any effects from the power of the high frequency envelope) (Figure 1A),

y=AH∥AH∥2  and s=[cos⁢(PL)⁢sin⁢(PL)](12)

FIGURE 1

Figure 1. (A) Preprocessing steps and (B) LPV-AR identification for CFC quantification. (A) The signal of interest is first band-pass filtered at the frequency bands L and H. After applying the Hilbert transform, the instantaneous amplitude and phase signals are used to estimate the LPV-AR based CFC. (B) LPV-AR model order selection procedure is realized as described in section “Linear parameter varying autoregressive model order selection.” LPV-AR identification is achieved by estimating the model coefficients. A new LPV-AR model is identified by rearranging the regressor matrix. CFC is quantified based on the residuals of both models.

In Eq. 12 we used the cosine and the sine of the phase of the L band to avoid the so-called “null phases” (Penny et al., 2008). Since the scheduling variable s is now a two-dimensional signal (i.e., s ∈ RN×2 compared to the one-dimensional case in Eq. 2), Eq. 2 becomes,

s=[cos⁢(PL)⁢sin⁢(PL)]=[s1⁢s2](13)

y⁢(n)=∑k=1pak⁢(sn)⁢y⁢(n-k)+ε⁢(n)(14)

Based on the polynomial basis, the AR model coefficients are now expressed as,

ak(sn)=θkTΨk=θk0+∑0≤i+j≤qθki⁢js1⁢nis2⁢nj,k=1:p(15)

where for notational simplicity s1n = s1(n) and s2n = s2(n). To make models with different polynomial orders more easily identifiable, we retained only the q-th order interactions. The final model representation can be expressed as in Eq. 8.

For illustration purposes, assume that p = 2, q = 2 then Eq. 14 can be written as,

y(n)=(θ10+θ102s2⁢n2+θ111s1⁢ns2⁢n+θ120s1⁢n2)y(n-1)+(θ20

+θ202s2⁢n2+θ211s1⁢ns2⁢n+θ220s1⁢n2)y(n-2)+ε(n)(16)

To estimate the CFC index between s and y we first estimate the coefficients based on the regularized solution of Eq. 10. The predicted variable of interest y^ and the residuals e of the LPV-AR model are computed as,

y^=Φ⁢wRLS,e=y-y^(17)

Then, we permute columnwise the order of s to destroy any dependencies between s and y and we recalculate the matrix Φ of Eq. 7, generating this way a new regressor matrix Φ0. The idea behind this procedure is that if the instantaneous phase s modulates the instantaneous amplitude y, rearranging its values will have a negative impact on the prediction accuracy of the LPV-AR model. On the other hand, if the coupling between the two is not significant, the predictive performance of the model will not be affected strongly. Based on these changes in predictive performance, we can quantify the CFC between s and y. Using the rearranged regressor matrix Φ0, the predicted signal of interest and the residuals of the LPV-AR model are given as,

y^0=Φ0⁢wRLS,e0=y-y^0,(18)

The CFC index between s and y is defined as (Figure 1B),

C⁢F⁢C=|log⁡(eT⁢ee0T⁢e0)|(19)

The higher the CFC index, the stronger the coupling between the instantaneous phase and amplitude of the signals of interest.

Linear parameter varying autoregressive model order selection

To estimate the CFC index, an optimal model structure should be first selected. The LPV-AR model complexity is essentially defined by the model order p and the polynomial order q. The higher the p and q, the better the fit to the observed data but also the lower the predictive performance and thus the generalization capabilities to unseen data. Model order selection revolves around the optimal selection of these two hyperparameters. A common practice is to use cross-validation (CV) or model order selection criteria like the Akaike Information Criterion (Akaike, 1974) or the Bayesian Information Criterion (Schwarz, 1978). For the purposes of CFC estimation, we developed the two-step procedure described below:

1. For a fixed and relatively large regularization parameter λ, we estimate the mean squared error (MSE) on the whole dataset for an ascending order of p and q values. The p and q values that achieve the lowest MSE are chosen as optimal.

2. In this step, the p and q values from step 1 are kept fixed and we optimize the regularization parameter. To determine a suitable value, we use the U-curve technique (Krawczyk-StanDo and Rudnicki, 2007), whereby the sum of the inverse of the norm of the regularized solution and the inverse of the residual norm is plotted in a log-log scale for different regularization parameters (i.e., λ = ). The regularization parameter that corresponds to the minimum of the U-curve is then selected as optimal, since it provides a good trade-off between the size of the regularized solution, and its fit to the given data.

The purpose of step 1 is mainly to obtain an appropriate value for the polynomial order q, since we found that, based on the simulations and real data described in the following sections, it affects the results. We also observed that all p values return a MSE for the same order q (see section “Simulations” – Figure 3). The selection of p is important but less critical at this step. In step 2, we fine tune the regularization parameter based on the p and q values acquired at step 1. This way, we make sure that, despite the selection of maybe a fairly complex model in terms of p at step 1, a suitable regularization parameter is applied. Apart from controlling model complexity, regularization ensures stability [i.e., poles inside the unit circle (Ljung, 1998)]. Narrowband signals are known to induce temporal instabilities on the AR models, because the roots of the signal generating system are located very close to the unit circle (Hall et al., 1983; Kostoglou and Lunglmayr, 2020). By applying regularization, the chances of obtaining unstable estimates are low. In addition, it provides more consistent model behavior especially when consecutive time samples are highly correlated. This could happen, for example, in lower frequencies that are oversampled (i.e., the bandpass-filtered lower frequency components are not down-sampled. LPV-AR estimation is always conducted in our case using the initial sampling rate).

FIGURE 2

Figure 2. Representative realizations of synthetic signals (Eq. 31) from (A) simulation set I (monophasic coupling) and (B) simulation set II (biphasic coupling). All depicted signals exhibit CFC between fL = 4 Hz and fH = 60 Hz. The sampling rate was set to Fs = 240 Hz and the signal to noise ratio was 40 dB.

FIGURE 3

Figure 3. Two step model order selection procedure for one representative realization of Eq. 31 and simulation set I. (A) The top panel depicts the MSE obtained during step 1 for different p and q values using a regularization parameter of λ = 10 as an initial value. The bottom panel depicts the U-curve from step 2, using the model order that achieved the smallest MSE [i.e., (p,q) = (10,2)] at step 1 (i.e., top panel). The regularization parameter that corresponds to the minimum of the U-curve was selected as optimal (i.e., λ = 1). (B) Similar as (A), however, the initial regularization parameter at step 1 was set to λ = 0.01.

Commonly used cross frequency coupling estimation methods

Here, we provide a brief overview of some of the most used methods for CFC estimation. These methods will be later compared with our proposed LPV-AR approach using both simulated and real EEG data.

Mean vector length

To assess CFC Canolty et al. (2006) defined a complex time-series with magnitude equal to the instantaneous high frequency amplitude and angle equal to the instantaneous low frequency phase,

z=AH⁢ej⁢PL(20)

The MVL is the absolute value of the mean of z,

M⁢V⁢L=||1N⁢∑n=1Nz⁢(n)||(21)

In the absence of CFC, vector points of z(n) in the complex plane exhibit uniform radial symmetry around zero and therefore the MVL takes values close to zero. If there is modulation, then this symmetry is lost and MVL obtains larger values.

Phase-locking value

The PLV (Lachaux et al., 1999; Cohen, 2008) is inspired by the well-known phase-locking concept frequently used in neuroscience. For the calculation of the PLV, the high frequency amplitude signal is Hilbert transformed and the phase of the corresponding analytic signal is extracted,

PAH⁢∡⁢H⁢i⁢l⁢b⁢e⁢r⁢t⁢(22)

PLV is then defined as,

P⁢L⁢V=||1N⁢∑n=1Ne⁢x⁢p⁢||(23)

When the two phase series PL and PAH are locked, then PLV is 1. In case of desynchronization (no CFC) PLV is 0.

Kullback–Liebler based modulation index

The Kullback–Liebler based MI by Tort et al. (2010) measures the deviation of the distribution of the high frequency amplitude, with respect to the phase of the low frequency, from the uniform distribution. Specifically, the low frequency phase is binned and the mean of the high frequency amplitude over each phase bin is computed. Each bin value is normalized by the sum over all bins,

P⁢(i)=<AH>(i)∑b⁢i⁢n=1N⁢b⁢i⁢n⁢s<AH>(b⁢i⁢n) for i=1⁢…⁢N⁢b⁢i⁢n⁢s(24)

where < AH > (i) refers to the mean high frequency amplitude that corresponds to the i-th phase bin. The default value for the number of bins is usually set to Nbins = 18. The Shannon entropy of P is given as,

H⁢(P)=-∑i=1N⁢b⁢i⁢n⁢sP⁢(i)⁢log⁡(25)

and the Kullback–Liebler distance (that measures the discrepancy between two distributions) is defined as,

DK⁢L=log⁡(N⁢b⁢i⁢n⁢s)-H⁢(P)(26)

where log(Nbins) is the maximum possible entropy value if the distribution of P is uniform. The MI is computed as,

M⁢I=DK⁢L/log⁡(N⁢b⁢i⁢n⁢s)(27)

If the mean amplitude is uniformly distributed over phase (i.e., no CFC), then MI is zero. As the distribution of P gets further away from the uniform distribution the Kullback–Liebler distance and the MI increases.

Generalized linear model

Penny et al. (2008) introduced the GLM index computed based on a multiple regression model, whereby the high frequency amplitude is expressed a linear combination of the cosine and the sine of the low frequency phase,

AH=[cos⁢(PL)⁢sin⁢(PL)⁢   1]⋅wG⁢L⁢M+ε(28)

where ε is additive Gaussian noise and 1 is a column of ones. wGLM can be estimated using the Least Squares solution. The proportion of variance explained by the model is used as the GLM index.

Simulations Signal generation

To validate our approach, we simulated two sets of 100 6 s-long signals that exhibited CFC between fL = 4 Hz and fH = 60 Hz at a sampling rate of Fs = 240 Hz. The exact steps followed are summarized below:

1. We bandpass-filtered white Gaussian noise at a center frequency of fL = 4 Hz and with a bandwidth of 1 Hz using a two-way least squares FIR filter as described in section “Linear parameter varying autoregressive based cross frequency coupling.” The obtained signal xL was then normalized to unit variance.

2. The amplitude of a sinusoid at fH = 60 Hz was modulated by the signal xL (of step 1) as follows (Figure 2),

SimulationsetI:   xH(n)=11+e-3⁢xL⁢(n)sin(2πfHFsn)(29)

SimulationsetII:   xH(n)=(11+e-3⁢xL⁢(n)-0.5)

sin⁡(2⁢π⁢fHFs⁢n)(30)

3. The final signal was computed as,

x⁢(n)=xH⁢(n)+xL⁢(n)+u⁢(n)(31)

where u is zero-mean white Gaussian noise with standard deviation σu = 0.5 (around 40 dB signal to noise ratio).

In the two simulation sets the amplitude of the high frequency is modulated differently by the phase of the low frequency. Specifically, in set I, the high frequency amplitude is maximum during the peaks of the low frequency, known as monophasic coupling, whereas in set II, it attains its maximum during the peaks but also the troughs of the low frequency component, known as biphasic coupling. For each simulation set, we generated different realizations of Eq. 31 by repeating steps 1–3 100 times (i.e., xL and u were produced based on different Gaussian white noise realizations). Some representative examples for both simulation sets can be found in Figure 2.

Quantifying cross frequency coupling using comodulograms

Comodulograms are coupling maps that depict the CFC strength as a function of the phase frequency and the amplitude frequency. For each realization of simulation sets I and II, we created comodulograms for amplitude frequencies ranging between 20 and 80 Hz (with a 2 Hz step and a bandwidth of 2 Hz, i.e., 20–22, 22–24 Hz, etc.) and phase frequencies between 2 and 10 Hz (with a 1 Hz step and a bandwidth of 0.5 Hz). LPV-AR based CFC was estimated at each amplitude and phase frequency band as described in sections “Linear parameter varying autoregressive based cross frequency coupling” and “Linear parameter varying autoregressive model order selection.” Model order selection was realized at each bin. To compare our approach with other methods, we analyzed the average comodulograms, over all 100 realizations, extracted in each case.

EEG data Experimental protocol and preprocessing

In total, 61-channel EEG signals were obtained from 8, originally right-handed, SCI participants (median ± SD age of 54 ± 18) with a neurological level of injury ranging from C1 to C7 as described in Ofner et al. (2019).1 A fixation cross along with a beep sound was presented on a black screen to signal the beginning of the trial. Class cues were displayed 2 s after the trial initialization and each trial lasted for 5 s. The participants were asked to attempt unilaterally arm/hand movements such as hand open, supination, pronation, palmar grasp and lateral grasp. Each attempted movement class consisted of 72 trials.

The recorded EEG signals were pre-processed using EEGLAB and Matlab. Bandpass filtering between 0.3 and 70 Hz was applied to the signals using a fourth order, zero-phase, Butterworth filter. Trials with impulsive noise were rejected using thresholding. Techniques based on abnormal joint probabilities and kurtosis (Ofner et al., 2019) were also employed in order to eliminate noisy trials. Blinks, saccades, and muscle movements were removed using independent component analysis.

To reduce the dimensionality of the data we defined the following nine regions of interest (ROI) in the sensor space: FZ (average of channels FCz and Fz), FL (average of channels F3, F1, FFC5h, FFC3h, FFC1h, FC5, FC3, FC1, FCC5h, FCC3h, and FCC1h), FR (average of channels F2, F4, FFC2h, FFC4h, FFC5h, FC2, FC4, FC6, FCC2h, FCC4h, and FCC6h), CZ (average of channels Cz and CPz), CL (average of channels C5, C3, C1, CCP5h, CCP3h, CCP1h, CP5, CP3, CP1, CPP5h, CPP3h, and CPP1h), CR (average of channels C2, C4, C6, CCP2h, CCP4h, CCP6h, CP2, CP4, CP6, CPP2h, CPP4h, and CPP6h), PZ (average of channels Pz and POz), PL (average of channels P5, P3, P1, and PPO1h), and PR (average of channels P2, P4, P6, and PPO2h). Note that in order to estimate CFC we first filtered the signals in the frequency bands of interest and then we averaged over different ROIs.

Cross frequency coupling estimation

For each subject and each attempted movement, we estimated CFC on a ROI-by-ROI and trial-by-trial basis and in the following frequency bands of interest: delta (D: [0.5 3] Hz), theta (T: [4 8]), alpha (A: [8.5 12]), beta (B: [12.5 30]), gamma (G: [30.5 60]), and high gamma (HG: [60 70]). The pair of phase and amplitude frequency bands [L,H] (i.e., CFCLH) were [D,T],[D,A],[D,B], [D,G],[D,HG],[T,A],[T,B],[T,G],[T,HG],[A,B], [A,G],[A,HG],[B,G],[B,HG],[G,HG]. Instead of using the whole range of the phase frequency band, we defined as the low frequency phase driver the sub-band with the highest power (by bandpass-filtering the signals in incremental steps of 0.5 and 1 Hz bandwidth and estimating their corresponding power). Therefore, for each trial we extracted 135 CFC indices (15 indices × 9 ROI). Trial-based CFC was computed using the same techniques investigated in the simulations section, along with the LPV-AR approach. LPV-AR model order selection and estimation was conducted for each trial and for each pair of low and high frequency bands. The extracted CFC indices were then used as features for movement attempt classification. We considered only data segments corresponding to the period between the cue onset and the end of the trial (3 s long).

Movement attempt classification using cross frequency coupling

We classified movement attempts using a multiclass shrinkage linear discriminant analysis classifier (Peck and Van Ness, 1982; Ofner et al., 2019). We focused on five types of movements: hand open, pronation, supination, palmar grasp and lateral grasp. We trained different classifiers for each subject and for each CFC estimation method. Trial-based classification performance was assessed using accuracy.

To reduce the number of features and improve classification performance, for each method and each subject, we applied a backward elimination feature selection scheme. Initially we trained a classifier with 135 CFC indices. To assess feature importance, we randomly shuffled the values of each feature in the testing sets and computed the average change in accuracy. The idea is that if a feature is important, rearranging its values will lead to a drop in accuracy performance. If the feature, however, is non-informative then the testing set accuracy will not be affected significantly. At each iteration, the feature with the lowest importance was removed and new classifiers were trained. This procedure was repeated until no features were left. Based on this method, the accuracy increased with the number of features, reached a maximum point, and then decreased again. The number of features that corresponded to this maximum point, was selected as optimal.

To validate the classification results and the feature selection procedure, we used both external and internal CV schemes. The external CV was used to quantify the generalization capabilities of the classifier on unseen data (i.e., 10% hold-out). The internal CV (fivefold) was used for feature selection, i.e., different classifiers were trained on each training set and feature importance was assessed based on the mean CV testing performance. To test for statistical significance and eliminate the hypothesis of possible overfitting, the class labels were permuted randomly and the whole feature selection and classification procedure was repeated using the same train/test/validation sets applied on the non-permuted data. To assess the overall performance in different hold out sets during the external CV, the whole internal/external CV scheme was employed 50 times by randomly resampling 10% of the data for hold-out. To predict the hold-out set, for each subset of features, the classifier was retrained based on all the data used for internal CV.

Comparison of different cross frequency coupling estimation methods based on the classification results

Our main assumption was that the CFC method that produced the highest classification performance would also provide more accurate CFC estimates. For a fair comparison, we investigated classification performances based on different number of features. This is also one of the reasons that we applied the feature selection procedure of section “Movement attempt classification using cross frequency coupling” for each CFC estimation method.

Results Simulations

In Figure 3 we illustrate the results obtained from the two-step model order selection procedure for a representative realization of simulation set I. The top panels depict the MSE obtained during step 1 for different p and q values using a fixed regularization parameter (i.e., λ = 10 and λ = 0.01 for Figures 3A,B, respectively). The bottom panels depict the U-curves extracted during step 2, using the model order that achieved the smallest MSE at step 1 (i.e., top panels). The regularization parameter that corresponds to the minimum of the U-curve was selected as optimal. The difference between Figures 3A,B is the value of the regularization parameter that was used for step 1. In both cases, however, we ended up with the same polynomial order of q = 2, which produces the smallest MSE for all p values. Therefore, step 1 is important in terms of selecting and appropriate q value, and step 2 is applied to regularize the complexity imposed by the selected p from step 1.

In Figure 4 we illustrate the average comodulograms obtained over all 100 realizations for both simulation sets and for different CFC estimation methods. The 2D comodulograms were averaged either over all phase frequencies or over all amplitude frequencies producing the results illustrated in Figure 5. For the simulation set I (i.e., the monophasic case), all methods were able to detect CFC between fL = 4 Hz and fH = 60 Hz. However, the LPV-AR methodology exhibited higher resolution detection (i.e., higher CFC at fL = 4 Hz and fH = 60 Hz and lower sidebands around these frequencies) and lower bias in frequencies where no coupling existed. For the simulation set II, the biphasic nature of the coupling led to distorted measures of CFC in most algorithms. The LPV-AR methodology achieved the best CFC detection results, followed by the MI technique. The GLM, PLV, and MVL were not able to detect the coupled phase and amplitude frequencies correctly. In terms of computational complexity, the LPV-AR model requires a larger number of operations and therefore its runtime is expected to be higher. The mean runtime was 0.26 ± 0.06 ms for MI, 0.31 ± 1.70 ms for GLM, 1.13 ± 0.48 ms for PLV, 0.06 ± 0.01 ms for MVL, and 84.52 ± 14.88 ms for LPV-AR.

FIGURE 4

Figure 4. Average comodulograms over 100 realizations for (A) simulation set I and (B) simulation set II. Each row corresponds to a different method (i.e., MI first row, GLM second row, PLV third row, MVL fourth row, and LPV-AR fifth row). For comparison purposes, the average comodulogram of each method was normalized by its maximum value over all phase and amplitude bins.

FIGURE 5

Figure 5. The average of the comodulograms depicted in Figure 4 for (A) simulation set I and (B) simulation set II along the amplitude (top panels) and the phase frequencies (bottom panel) for different CFC methods (normalized to their maximum value). The LPV-AR technique exhibits higher resolution in detecting CFC (i.e., higher CFC at fL = 4 Hz and fH = 60 Hz and lower sidebands around these frequencies) and lower bias in frequencies where no coupling exists.

EEG data

The average accuracy plots (over all participants and over all hold out sets) for the five-class EEG classification problem (i.e., hand open vs. palmar grasp vs. lateral grasp vs. pronation vs. supination) as a function of the number of features can be found in Figure 6. The LPV-AR approach achieved the highest accuracy (Internal CV: 58 ± 5% – External CV: 54 ± 5%), followed by the MI (Internal CV: 55 ± 5% – External CV: 51 ± 5%), the GLM (Internal CV: 53 ± 4% – External CV: 49 ± 5%), the PLV (Internal CV: 50 ± 5% – External CV: 47 ± 6%), and finally the MVL (Internal CV: 45 ± 3% – External CV: 40 ± 3%). Based on the Kruskal–Wallis test, for 47 features (number of features that correspond to the maximum internal CV accuracy for the LPV-AR model) significant differences in accuracy performance were detected between methods. By applying the Dunn’s post hoc test for multiple comparisons analysis, the LPV-AR model was found to perform significantly better than the other methods in both external and internal CV (p < 0.05).

FIGURE 6

Figure 6. Average accuracy over all participants for the five-class classification problem (i.e., hand open vs. palmar grasp vs. lateral grasp vs. pronation vs. supination) as a function of the number of features during (A) the internal cross-validation and (B) the external cross-validation step. Different line colors correspond to different CFC estimation methods. Dashed lines represent accuracy levels obtained by training the classifier on data with permuted labels. The theoretical chance accuracy level is 20%. In (A), the maximum classification accuracy for each CFC method is indicated with a black circle.

Cross frequency coupling feature