Re-thinking EEG-based non-invasive brain interfaces: modeling and analysis
Gaurav Gupta Sergio Pequito Paul Bogdan












































































Re-thinking EEG-based non-invasive brain interfaces: modeling and analysis
    Gaurav Gupta1,
Sérgio Pequito2 and
Paul Bogdan1
1Ming Hsieh Department of Electrical Engineering, University of Southern California, Los Angeles, CA, USA

Email: {ggaurav, pbogdan}@usc.edu
2Department of Industrial and Systems Engineering, Rensselaer Polytechnic Institute, Troy, NY, USA

Email: goncas@rpi.edu

    December 30, 2023
===========================================================================================================================================================================================================================================================================================================================





Brain interfaces are cyber-physical systems that aim to harvest information from the (physical) brain through sensing mechanisms, extract information about the underlying processes, and decide/actuate accordingly. Nonetheless, the brain interfaces are still in their infancy, but reaching to their maturity quickly as several initiatives are released to push forward their development (e.g., NeuraLink by Elon Musk and `typing-by-brain' by Facebook). This has motivated us to revisit the design of EEG-based non-invasive brain interfaces. Specifically, current methodologies entail a highly skilled neuro-functional approach and evidence-based a priori knowledge about specific signal features and their interpretation from a neuro-physiological point of view. Hereafter, we propose to demystify such approaches, as we propose to leverage new time-varying complex network models that equip us with a fractal dynamical characterization of the underlying processes. Subsequently, the parameters of the proposed complex network models can be explained from a system's perspective, and, consecutively, used for classification using machine learning algorithms and/or actuation laws determined using control system's theory. Besides, the proposed system identification methods and techniques have computational complexities comparable with those currently used in EEG-based brain interfaces, which enable comparable online performances. Furthermore, we foresee that the proposed models and approaches are also valid using other invasive and non-invasive technologies. Finally, we illustrate and experimentally evaluate this approach on real EEG-datasets to assess and validate the proposed methodology. The classification accuracies are high even on having less number of training samples.




brain interfaces; spatiotemporal; fractional dynamics; unknown inputs; classification; motor prediction























§ INTRODUCTION



We have recently testimony a renewed interest in invasive and non-invasive brain interfaces. Elon Musk has released the NeuraLink initiative <cit.> that aims to develop devices and mechanisms to interact with the brain in a symbiotic fashion, thus merging the artificial intelligence with the human brain. The potential is enormous since it would ideally present a leap in our understanding of the brain, and an unseen enhancement of its functionality. Alternatively, Facebook just announced the `Typing-by-Brain' project <cit.> that gathered a team of 60 researchers whose target is to be capable of writing 100 words per minute that contrasts with the state-of-the-art of 0.3 to 0.82 words per minute assuming an average of 5 letters per word. Towards this goal, Facebook has invested in developing new non-invasive optical imaging technology that is five times faster and portable with respect to the one available on the market and would possess increased spatial and temporal resolution. Nonetheless, these are just some of the initiatives among others by some big Silicon Valley players that want to commercialize brain technologies <cit.>.

Part of the motivation for the `hype' in the use of neurotechnologies – both invasive and non-invasive brain interfaces – is due to their promising application domains <cit.>: (i) replace, i.e., the interaction of the brain with a wheelchair or a prosthetic device to replace a permanent functionality loss, (ii) restore, i.e.,  to bring to its normal use some reversible loss of functionality such as walking after a severe car accident or limb movement after a stroke, (iii) enhance, i.e.,  to enable to outperform in a specific function or task, as for instance an alert system to drive for long periods of time while attention is up; and (iv) supplement as in equipping one with extra functionality, as a third arm to be used during surgery. Notwithstanding, these are just some of the (known) potential uses of neurotechnology.

Despite the developments and promise of future applications of brain interfaces (some of which we cannot currently conceive), we believe that current approaches to both invasive and non-invasive brain interfaces can greatly benefit from cyber-physical systems (CPS) oriented approaches and tools to increase their efficacy and resilience. Hereafter, we propose to focus on non-invasive technology relying on electroencephalogram (EEG) and revisit it through a CPS lens. Moreover, we believe that the proposed methodology can be easily applicable to other technologies, e.g., electromagnetic fields (magnetoencephalography (MEG) <cit.>, and the  hemodynamic responses associated to neural activity, e.g. functional magnetic resonance imaging (fMRI) <cit.>, and functional
near-infrared spectroscopy (fNIRS) <cit.>). Nonetheless, these technologies present several drawbacks compared to EEG, e.g., cost, scalability, and user comfort, which lead us to focus on EEG-based technologies. Similar argument can be applied in the context of non-invasive versus invasive technologies that require patient surgery.



Traditional approach to EEG neuroadaptive technology consists of  proceeding through the following steps <cit.>: (a) signal acquisition (in our case, measurement of the EEG time-series); (b) signal processing (e.g., filtering with respect to known error sources); (c) feature extraction (i.e., an artificial construction of the signal that aims to capture quantities of interest); (d) feature translation (i.e., classification of the signal according to some plausible neurophysiological hypothesis); and (e) decision making (e.g., provide instructions to the computer to move a cursor, or a wheelchair to move forward) – see also Figure <ref> for an overview diagram. 

In this paper, we propose to merge steps (b) and (c) motivated by the fact that these often accounts for spatial or temporal properties, and are only artificially combined in a later phase of the pipeline, i.e., at step (d) of feature translation. First, we construct a time-varying complex network where the activity of nodes (or EEG signal) have long-range memory and the edges accounts for inter-node spatial dependence. Thus, we argue that the previous approach discards several spatial-temporal properties that can be weighted for signal processing and feature extraction phases. In other words, current EEG brain interfaces require one to have an understanding of the different regions of the brain responsible. For instance, regions for motor or visual actions, as well as artificial frequency bands that are believed to be more significant for a specific action (also known as evidence-based). Besides, one needs to understand and anticipate the most likely causes noise/artifacts in the EEG data collected and filter out entire frequency bands, which possibly compromises phenomena of interest not being available for post-processing. Instead, we propose a modeling capability to enable the modeling of long-range memory time-series that at the same time accounts for unknown stimuli, e.g., artifacts or inter-region communication.





 §.§ Related Work and Novel Contributions


We put forward that the ability to properly model the EEG time-series with complex network models, that account for realistic setups, can enable the brain interfaces methods to get transformed from detectors to decoders. In other words, we do not want to solely look for the existence of a peak of activity in a given band that is believed to be associated with a specific action. But we want to decompose the bunch of signals into different features, i.e., parameters of our complex network model, that are interpretable. Thus, allowing us to understand how different regions communicate during a specific action/task by representing them as nodes of the complex network and estimating the dependence via coupling matrix. The node activities are assumed to be evolving at different time-scales and in the general setup of presence of external stimuli which could either be unwanted noise or external process driving the system. In engineering, this will enable us to depart from a skill dependent situation to general context analysis, which will enhance the resilience of the approaches for practical nonsurgical brain interfaces. Besides, it will equip bioengineers, neuroscientists, and physicians with an exploratory tool to pursue new technologies for neuro-related diagnostics and treatments, as well as neuro-enhancement.































The proposed approach departs from those available in the literature, see <cit.> and references therein. In fact, to the best of authors' knowledge, in the context of noninvasive EEG-based technologies, <cit.> is the only existing work that explores fractional-order models in the context of single-channel analysis, which contrasts with the spatiotemporal modeled leveraged in this paper that copes with multi-channel analysis. Furthermore, the methodology presented in this paper also accommodate unknown stimuli <cit.>. For which efficient algorithms are proposed and leveraged hereafter to simultaneously retrieve the best model that conforms with unknown stimuli, and separating the unknown stimuli from the time-series associated with brain activity.  Our methods are as computationally efficient and stable as least-squares and spectral analysis methods used in a spatial and temporal analysis, respectively; thus, suitable for online implementation in nonsurgical brain interfaces.








The main contributions of the present paper are those of leveraging some of the recently proposed methods to develop new modeling capabilities for the EEG based neuro-wearables that are capable of enhancing the signal quality and decision-making. Furthermore, the parametric description of these models provides us with new features that are biologically motivated and easier to translate in the context of brain function associated with a signal characterization, and free of signal artifacts. Thus, making the brain-related activity interpretable, which leads to resilient and functional nonsurgical brain interfaces.



 §.§ Paper Organization


The remaining of the paper is organized as follows. Section <ref> introduces the complex network model considered in this paper and the main problem studied in this manuscript. Also we will see the description of the employed method for feature selection and then classification techniques. In Section <ref>, we present an elaborated study on the different datasets taken from the BCI competition <cit.>.




§ RE-THINKING EEG-BASED NON-INVASIVE 
 BRAIN INTERFACES



Brain interfaces aim to address the following problem.

Is it possible to classify a specific cognitive state, e.g., motor task or its imagination, by using measurements collected with a specific sensing technology that harvest information about brain activity?

In the current manuscript, we revisit this problem in the context of brain-computer interfaces (BCI), when dealing with EEG-based noninvasive brain interfaces. Towards this goal, we review the currently adopted procedures for solving this problem (see Figure <ref> for an overview), and propose a systems' perspective that enables to enhance the BCI reliability and resilience. Therefore, in Section <ref> we provide a brief overview of the EEG-based technology and the connection with the brain-areas' function associated with studies conducted in the past. Next, in Section <ref>, we introduce the spatiotemporal fractional model under unknown stimuli. This will be the core of the proposed approach in this manuscript to retrieve new features for classification, and, subsequently, enhancing brain interfaces capabilities. In Section <ref>, we describe how to determine the system model's parameters, and in Section <ref> we describe how to interpret them for the classification task.





 §.§ EEG-based Technology for Brain Interfaces: a brief overview



EEG enables the electrophysiological monitoring of space-averaged synaptic source activity from millions of neurons occurring at the neocortex level. The EEG signals have a poor spatial resolution but high temporal resolution, since the electrical activity generated at the ensemble of neurons level arrives at the recording sites within milliseconds. The electrodes (i.e., sensor) are placed over an area of the brain of interest, being the most common the visual, motor, sensory, and pre-frontal cortices. Usually, they follow standard montages – the International 10-20 system is depicted in Figure <ref>.

Most of the activity captured by the EEG electrodes is due to the interactions between inhibitory interneurons and excitatory pyramidal cells, which produces rhythmic fluctuations commonly referred to as oscilations. The mechanisms that generate those oscillations is not yet completely understood, but it has been already identified that some `natural oscillations' provide evidence of activity being `processed' in certain regions of the brain at certain `frequencies'. Therefore, oscillatory behavior of human brain is often partitioned in bands (covering a wide range of frequencies decaying as 1/f in power): (i) δ-band (0.5-3Hz); (ii) θ-band (3.5-7Hz); (iii) α-band (8-13Hz); (iv) β-band (14-30Hz); and (v) γ-band (30-70Hz). Furthermore, there has been some evidence that activity in certain bands is associated with sensory registration, perception, movement and cognitive processes related to attention, learning and memory <cit.>. Notwithstanding, such associations are often made using correlation and/or coherence techniques that only capture relationships between specific channels. But such methods are not able to assess the causality between signals that enables forecasting on the signal evolution, captured by a model-based representation as we propose to do hereafter. 

Different changes in the signals across different bands are also used to interpret the event-related potentials (ERPs) in the EEG signals, i.e., variations due to specific events – see <cit.> for detailed analysis. In the context of sensory-motor data used in the current manuscript to validate the proposed methodology, sensorimotor rhythms (SMRs) are often considered. These represent oscillations that are recorded over the posterior frontal and anterior parietal areas of the brain, i.e., over the sensorimotor cortices (see Figure <ref>). SMRs occur mainly in the α-band (for sensors located on the top of the motor cortices), and on beta and lower gamma for those on the sensorimotor cortices <cit.>. Consequently, these have been used as a default feature for classification of motor-related execution and only the imagination of performing a motor task. Notwithstanding, the spatiotemporal modeling is simultaneously captured through direct state-space modeling that enables the system's understanding of the dynamics of the underlying process. In addition, it provides a new set of attributes that can be used to improve feature translation, i.e., classification. 




 §.§ Spatiotemporal Fractional Model with Unknown Stimuli




A multitude of complex systems exhibit long-range (non-local) properties, interactions and/or dependencies (e.g., power-law decays in memories). Example of such systems includes Hamiltonian systems, where memory is the result of stickiness of trajectories in time to the islands of regular motion <cit.>. Alternatively, it has been rigorously confirmed that viscoelastic properties are typical for a wide variety of biological entities like stem cells, liver, pancreas, heart valve, brain, muscles <cit.>, suggesting that memories of these systems obey the power law distributions. These dynamical systems can be characterized by the well-established mathematical theory of  fractional calculus <cit.>, and the corresponding systems could be described by fractional differential equations <cit.>.  However, it is until recently that fractional order system (FOS) starts to find its strong position in a wide spectrum of applications in different domains. This is due to the availability of computing and data acquisition methods to evaluate its efficacy in terms of capturing the underlying system states evolution. 

Subsequently, we consider a linear discrete time fractional-order dynamical model under unknown stimuli (i.e., inputs) described as follows:

    Δ^αx[k+1]    =    Ax[k] + Bu[k] 
    
    y[k]    =    Cx[k],

where x∈ℝ^n is the state, u ∈ℝ^p is the unknown input and y ∈ℝ^n is the output vector. The differencing operator Δ is used as the discrete version of the derivative, for example Δ^1x[k] = x[k] - x[k-1]. As evident, the difference order of 1 has only one-step memory, and hence the classic linear-time invariant models are not able to answer the long-range memory property of several physiological signals as discussed before. On the other hand, the expansion of fractional-order derivative  in the discretized setting <cit.> for any ith state (1≤ i≤ n) can be written as


    Δ^α_ix_i[k] = ∑_j=0^kψ(α_i,j) x_i[k-j],

where α_i is the fractional order corresponding to the ith state and ψ(α_i,j) = Γ(j-α_i)/Γ(-α_i)Γ(j+1) with Γ(.) denoting the gamma function. We can observe from (<ref>) that fractional-order derivate provide long-range memory by including all x_i[k-j] terms. A quick comparison between the prediction accuracy of fractional-order derivative model and linear-time invariant model is shown in Figure <ref>. The fractional-order model can cope with sudden changes in the signals while the linear model cannot.



We can also describe the system by its matrices tuple (α,A, B, C) of appropriate dimensions. The coupling matrix A represents the spatial coupling between the states across time while the input coupling matrix B determines how inputs are affecting the states. In what follows, we assume that the input size is always strictly less than the size of state vector, i.e., p < n. 

Having defined the system model, the system identification, i.e., estimation of model parameters, from the given data is an important step. It becomes nontrivial when we have unknown inputs since one has to be able to differentiate which part of the evolution of the system is due to its intrinsic dynamics and what is due to the unknown inputs. Subsequently, the analysis part that we need to address is that of system identification from the data, as described next.



 §.§ Data driven system identification



The problem consists of estimating α, A and inputs {u[k]}_t^t+T-2 from the given limited observations y[k], k = [t, t + T -1], which due to the dedicated nature of sensing mechanism is same as x[k] and under the assumption that the input matrix B is known. The realization of B can be application dependent and is computed separately using experimental data – as we explore later in the case study, see Section <ref>. For the simplicity of notation, let us denote z[k] = Δ^αx[k+1] with k chosen appropriately. The pre-factors in the summation in (<ref>) grows as ψ(α_i,j) ∼𝒪(j^-α_i-1) and, therefore, for the purpose of computational ease we would be limiting the summation in (<ref>) to J values, where J>0 is sufficiently large. Therefore, z_i[k] can be written as

    z_i[k] = ∑_j=0^J-1ψ(α_i,j)x[k+1-j],

with the assumption that x[k], u[k] = 0 for k≤ t-1. Using the above introduced notations and the model definition in (<ref>), the given observations can be written as

    z[k] = Ax[k] + Bu[k] + e[k],

where e ∼𝒩(0,Σ) is assumed to be  Gaussian noise independent across space and time. For simplicity we would assume that Σ = σ^2 I. Also, let us denote the system matrices as A = [a_1,a_2,,a_n]^T and B = [b_1,b_2,,b_n]^T. The vertical concatenated states and inputs during an arbitrary window of time as  X_[t-1,t+T-2] = [x[t-1],x[t],,x[t+T-2]]^T, U_[t-1,t+T-2] = [u[t-1],u[t],,u[t+T-2]]^T respectively, and for any ith state we have Z_i,[t-1,t+T-2] = [z_i[t-1],z_i[t],,z_i[t+T-2]]^T. For the sake of brevity we would be dropping the time horizon subscript from the above matrices as it is clear from the context.

Since the problem of joint estimation of the different parameters is highly nonlinear, we proceed as follows: (i) we estimate the fractional order α using the wavelet technique described in <cit.>; and (ii) with α known, the z in equation (<ref>) is computed under the additional assumption that the system matrix B is known. Therefore, the problem now reduces to estimate A and the inputs {u[k]}_t^t+T-2. Towards this goal, we exploit the algorithm similar to expectation-maximization (EM) <cit.> from <cit.>. Briefly, the EM algorithm is used for maximum likelihood estimation (MLE) of parameters subject to hidden variables. Intuitively, in our case, in Algorithm 1, we estimate A in the presence of hidden variables or unknown unknowns {u[k]}_t^t+T-2. Therefore, the `E-step' is performed to average out the effects of unknown unknowns and obtain an estimate of u, where due to the diversity of solutions, we control the sparsity of the inputs using the parameter λ'. Subsequently, the `M-step' can then accomplish MLE estimation to obtain an estimate of A.

It was shown theoretically in <cit.> that the algorithm is convergent in the likelihood sense. It should also be noted that the EM algorithm can converge to saddle points as exemplified in <cit.>. The Algorithm <ref> being iterative is crucially dependent on the initial condition for the convergence. We will see in Section <ref> that the convergence is very fast making it suitable for online estimation of parameters.














 §.§ Feature Translation (Classification)



The unprocessed EEG signals coming from the sensors although carrying vital information may not be directly useful for making the predictions. However, by representing the signals in terms of parametric model (α, A) and the unknown signals as we did in the last section, we can gain better insights. The parameters of the model being representative of the original signal itself can be used to make a concise differentiation. 

The A matrix represents how strong is the particular signal and how much it is affecting/being affected by the other signals that are considered together. While   performing or imagining particular motor tasks, certain regions of the brain gets more activated than others. Simultaneously, the inter-region activity also changes. Therefore, the columns of A which represent the coefficients of the strength of a signal affecting other signals can be used as a feature for classification of motor tasks. In this work, we will be considering the machine learning based classification techniques like logistic regression and Support Vector Machines (SVM) <cit.>. The other classification techniques c The choice of kernels would vary from simple `linear' to radial basis function (RBF), i.e., k(x_i,x_j) = e^-γ(x_i-x_j)^2. The value of parameters of the classifier and possibly of the kernels are determined using the cross-validation. The range of parameters in the cross-validation are from 2^-5,,2^15 for γ and 2^-15,,2^3 for C = 1/ λ, both in the logarithmic scale, where λ is the regularization parameter which appears in optimization cost of the classifiers <cit.>.






§ CASE STUDY



We will now illustrate the usefulness of the  dynamic model with unknown inputs in the context of classification for BCI. We have considered two datasets from the BCI competition <cit.>. The datasets were selected on the priority of larger number of EEG channels and number of trials for training. The available data is split into the ratio of 60% and 40% for the purpose of training and testing, respectively.



 §.§ Dataset-I



We consider for validation the dataset labeled `dataset IVa' from BCI Competition-III <cit.>. The recording was made using BrainAmp amplifiers and a 128 channel electrode cap and out of which 118 channels were used. The signals were band-pass filtered between 0.05 and 200 Hz and then digitized at 1000 Hz. For the purpose of this study we have used the downsampled version at 100 Hz. The dataset for subject ID `al' is considered, and it contains 280 trials. The subject was provided a visual cue, and immediately after asked to imagine two motor tasks: (R) right hand, and (F) right foot.



  §.§.§ Sensor Selection and Modeling

To avoid the curse-of-dimensionality, instead of considering 118 sensors available, which implies the use of 118× 118 dynamics entries for classification, only a subset of  9 sensors is considered. Specifically, only the sensors indicated in Figure <ref> are selected on the basis that only hand and feet movements need to be predicted, and only a 9× 9 dynamics matrix and 9 fractional order coefficients are required for modeling the fractional order system. Besides, these sensors are selected because they are close to the region of the brain known to be associated with motor actions. 





  §.§.§ System Identification and Validation

The model parameters (α,A) and the unknown inputs are estimated by using the Algorithm <ref>. As mentioned before, the performance of the algorithm being iterative is dependent on the choice of the initial conditions. For the current case, we have observed that the algorithm converges very fast, and even a single iteration is enough. The convergence of mean squared error in the M-step of Algorithm <ref> for one sample from dataset-I is shown in Figure <ref>. This shows that the choice of initial conditions are fairly good. The one step and five step prediction of the estimated model is shown in Figure <ref>. It is evident that the predicted values for one step very closely follow the actual values. There are some differences between the actual and predicted values for five step prediction.





  §.§.§ Discussion of the results
 
The most popular features used in the motor-imagery based BCI classification relies on exploiting the spectral information. The features used are the band-power which quantifies the energy of the EEG signals in certain spectrum bands <cit.>. The motor cortex of the brain is known to be affecting the energy in the bands namely, α and β as discussed in Section <ref>. While it happens that unwanted signal energy is captured in these bands as well while performing the experiments, for example neck movement, other muscle activities etc. The filtering of these so called `unwanted' components from the original signal is a challenging task using the spectral techniques as they often share the same band. 

We used a different approach to deal with these unknown unknowns in Section <ref>. The magnitude spectrum of the original EEG signal and on removing the estimated unknown inputs is shown in Figure <ref>. It should be observed that the original signal and the signal upon removing the unknown inputs have significant energy in the α and β bands. The unknown inputs behave similar to the white noise which is evident from their Gaussian probability distribution (PDF) as shown in Figure <ref>. The inputs are not mean zero but their PDF is centered around a mean value of approximately 58.






The model parameters (α, A) are jointly estimated with the unknown inputs using Algorithm <ref>, therefore the effect of the inputs is inherently taken care of in the parameters. The structure of matrix A for two different labels is shown in Figure <ref>. We have used the sensors C_3 and C_1 which are indexed as 3 and 4, respectively in Figure <ref>. It is apparent from Figure <ref> that the columns corresponding to these sensors have different activity and hence deem to be fair candidates for the features to be used in classification. Therefore, the total number of features are 2× 9=18.






 §.§ Dataset-II



A 118 channel EEG data from BCI Competition-III, labeled as `dataset IVb' is taken <cit.>. The data acquisition technique and sampling frequencies are same as in dataset of the previous subsection. The total number of labeled trials are 210. The subjects upon provided visual cues were asked to imagine two motor tasks, namely (L) left hand and (F) right foot.



  §.§.§ Sensor Selection and Modeling

Due to the small number of training examples, we have again resorted to select the subset of sensors for the model estimation as we did for the dataset-I in the previous section. Since the motor tasks were left hand and feet, therefore we have selected the sensors in the right half of the brain and close to the region which is known to be associated with hand and feet movements as shown in Figure <ref>. We will see in the final part of this section that selecting sensors based on such analogy helps not only in reducing the number of features, but also to gain better and meaningful results.









  §.§.§ System Identification and Validation

After performing the estimation of the model (α, A) and the unknown inputs using the subset of sensors, we can see the similar performance of the model on dataset-II as was in dataset-I. The convergence of mean squared error in the M-step of Algorithm <ref> for one sample from dataset-II is shown in Figure <ref>. The one step and five step predictions are shown in Figure <ref>. The model prediction follows closely the original signal.








  §.§.§ Discussion of the results
 
The spectrum of the original EEG signal at channel CFC_2 and its version with unknown inputs removed are shown in Figure <ref>. The spectrum shows peaks in the α and β bands. We witness a similar observation as before that both of the signals share the same band and hence making it difficult to remove the effects of the unwanted inputs. The unknown inputs resembles that of white noise and the PDF is close to Gaussian distribution with mean centered at around 48.

The estimated A matrix from Algorithm <ref> is shown in Figure <ref> for two different labels. Out of all 13 sensors, the sensors CCP_2 and CCP_4 which are indexed as 10 and 11 in the matrix have striking different activity. The columns corresponding to these two sensors seem good choice for being the features for classification. Therefore, the total number of features are 2× 13 = 26 for this dataset. Next, we discuss the classification accuracy for both the datasets.







 §.§ Classification Performance


Finally, the performance of the classifiers using the features explained for both the datasets are shown in Figure <ref>. The classifiers are arranged in the order of complexity from left to right with logistic regression (lR) and linear kernel being simplest and SVM with RBF kernel being most complex. The performance plot parallels the classic machine learning divergence curve for both the datasets. The accuracy for training data increases when increasing the classification model complexity while it reduces for the testing data. This is intuitive because a complex classification model would try to better classify the training data. But the performance of the test data would reduce due to overfitting upon using the complex models. We have very few training examples to build the classifier and hence such trend is expected. The performance of the classifiers for both the datasets are fairly high which reflects the strength of the estimated features. We can see a 87.6% test accuracy for dataset-I and 85.7% for dataset-II. While these accuracies depend a lot on the cross-validation numbers and other factors like choice of classifier which can be better tuned to get higher numbers. 

For both the datasets we have seen that the proposed methodology efficiently extracts the features which serves as good candidate to differentiate the imagined motor movements. By implicitly removing the effects of the unwanted stimuli, the coefficients of the coupling matrix A are shown to be sufficient for discriminating relation between various EEG signals which are indicative of the motor movements. The testing accuracies are high which indicate the good quality of the extracted features.




§ CONCLUSION



We have revisited the EEG-based noninvasive brain interfaces feature extraction and translation from a cyber-physical systems' lens. Specifically, we leveraged spatiotemporal fractional-order models that cope with the unknown inputs. The fractional-order models provide us the dynamic coupling changes that rule the EEG data collected from the different EEG sensors, and the fractional-order exponents capture the long-term memory of the process. Subsequently, unknown stimuli is determined as the external input that least conforms with the fractional-order model. By doing so, we have filtered-out from the brain EEG signals the unknown inputs, that might be originated in the deeper brain structures. The presence of unknown stimuli is possibly the result of the structural connectivity of the brain that crisscrosses different regions, or due to artifacts originated in the muscles (e.g., eye blinking or head movement). As a consequence, the filtered signal does not need to annihilate an entire band in the frequency domain, thus keeping information about some frequency regions of the signal that would be otherwise lost. 

We have shown how the different features obtained from the proposed model can be used towards rethinking the EEG-based noninvasive interfaces. In particular, two datasets used in BCI competitions were used to validate the performance of the methodology introduced in this manuscript, which is compatible with some of the state-of-the-art performances while requiring a relatively small number of training points. We believe that the proposed methodology can be used within the context of different neurophysiological processes and corresponding sensing technologies. Future research will focus on leveraging additional information from the unknown inputs retrieved to anticipate specific artifacts and enable the deployment of neuro-wearables in the context of real-life scenarios.  Furthermore, the presented methodology can be used as an exploratory tool by neuroscientists and physicians, by testing input and output responses and tracking their impact in the unknown inputs retrieved by the algorithm proposed; in other words, one will be able to systematically identify the origin and dynamics of stimulus across space and time. Finally, it would be interesting to explore the proposed approach in the closed-loop context, where the present models would benefit from control-like strategies to enhance the brain towards certain tasks or attenuate side effects of certain neurodegenerative diseases or disorders.





§ ACKNOWLEDGMENT



The authors are thankful to the reviewers for their valuable feedback. G.G. and P.B. gratefully acknowledge the support by the U.S. Army Defense Advanced Research Projects Agency (DARPA) under grant no. W911NF-17-1-0076, DARPA Young Faculty Award under grant no. N66001-17-1-4044, and the National Science Foundation under CAREER Award CPS-1453860 support.








IEEEtran


