## 2 Methods
### Kernel-Based Relevance Analysis  
When characterizing EEG data through effective brain connectivity measures for BCI-related applications, two common and related issues can arise. First, all to all channel connectivity analyses result in a large number of features, many of which may not provide useful information to discriminate between the conditions of the BCI paradigm of interest [10]. This can add noise and complexity to any subsequent analysis stage. Second, estimating such a large number of pairwise channel connectivities can be computationally expensive, especially for measures such as TE and single-trial TE\({}^{\theta}\)[8], which can hinder their inclusion in practical BCI systems. Both problems could be addressed by identifying the set of pairwise channel connectivities that are relevant to discriminate between specific conditions, which would also lead to a clearer neurophysiological interpretation of the obtained results [6; 10]. To that end, we explore a relevance analysis strategy based on centered kernel alignment (CKA). CKA allows quantifying the similarity between two sample spaces by comparing two characterizing kernel functions [29]. First, assume we have a feature matrix \(\mathbf{\Phi}\ \in\ \mathbb{R}^{N\times P}\), and a corresponding vector of labels \(I\ \in\ \mathbb{Z}^{N}\), with \(N\) the number of observations and \(P\) the number of features. For the case of connectivity-based EEG analysis, each element in \(\mathbf{\Phi}\) holds a connectivity value for a pair of channels, with each row of \(\mathbf{\Phi}\) containing multiple connectivity values (features) estimated for a single trial or observation. The corresponding element in \(I\) holds a label identifying the condition associated to that trial. Next, we define two kernel matrices \(\mathbf{K_{\Phi}}\ \in\ \mathbb{R}^{N\times N}\) and \(\mathbf{K_{I}}\ \in\ \mathbb{R}^{N\times N}\). The first matrix holds elements \(k_{ij}^{\mathbf{\Phi}}=\kappa_{\mathbf{\Phi}}(\mathbf{\phi}_{i},\mathbf{\phi}_{j})\) with \(\mathbf{\phi}_{i},\ \mathbf{\phi}_{j}\ \in\ \mathbb{R}^{P}\) row vectors belonging to \(\mathbf{\Phi}\), and  
\[\kappa_{\mathbf{\Phi}}(\mathbf{\phi}_{i},\mathbf{\phi}_{j};\sigma)=\exp\left(-\frac{d^ {2}(\mathbf{\phi}_{i},\mathbf{\phi}_{j})}{2\sigma^{2}}\right) \tag{13}\]  
a radial basis function (RBF) kernel [37], where \(d^{2}(\cdot,\cdot)\) is a distance operator and \(\sigma\ \in\ \mathbb{R}^{+}\) is the kernel's bandwidth. The second matrix has elements \(k_{ij}^{I}=\kappa_{I}(l_{i},l_{j})\) with \(l_{i},\ l_{j}\ \in\ I\), and  
\[\kappa_{I}(l_{i},l_{j})=\delta(l_{i}-l_{j}), \tag{14}\]  
a dirac kernel, where \(\delta(\cdot)\) stands for the Dirac delta. Then, the CKA can be estimated as:  
\[\hat{\rho}(\tilde{\mathbf{K}}_{\mathbf{\Phi}},\tilde{\mathbf{K}}_{I})=\frac{\langle\tilde{ \mathbf{K}}_{\mathbf{\Phi}},\tilde{\mathbf{K}}_{I}\rangle_{\mathbb{F}}}{(\langle\tilde{ \mathbf{K}}_{\mathbf{\Phi}},\tilde{\mathbf{K}}_{\mathbf{\Phi}}\rangle_{\mathbb{F}}\langle\tilde{ \mathbf{K}}_{I},\tilde{\mathbf{K}}_{I}\rangle_{\mathbb{F}})^{1/2}}, \tag{15}\]  
In this work, we use NMMs to generate interacting time series with known oscillatory properties in order to test the performance of the proposed phase TE estimator. In particular, we test our proposal in terms of its ability to detect directed interactions for different levels of coupling strength, under the presence of noise and signal mixing, and for bidirectional narrowband couplings. We proceed as follows: first, we set the model parameters describing Areas 1 and 2 as in [40], so as to generate signals with power spectrums peaking in the \(\alpha\) (8 Hz-12 Hz) and lower \(\beta\) bands (12-20 Hz), as depicted in Figure 1C. Then, in order to generate unidirectionally coupled signals, with interactions from \(\mathbf{x}\) to \(\mathbf{y}\), we set the parameter \(C_{21}\) to 0 for all simulations. Also, the parameter \(\
u\) is set to 20 ms, and the extrinsic inputs \(\mathbf{p}_{1}\) and \(\mathbf{p}_{2}\) are modeled as Gaussian noise [40].