On the Assessment of Functional Connectivity in an Immersive Brain-Computer Interface During Motor Imagery

New trends on brain-computer interface (BCI) design are aiming to combine this technology with immersive virtual reality in order to provide a sense of realism to its users. In this study, we propose an experimental BCI to control an immersive telepresence system using motor imagery (MI). The system is immersive in the sense that the users can control the movement of a NAO humanoid robot in a first person perspective (1PP), i.e., as if the movement of the robot was his/her own. We analyze functional brain connectivity between 1PP and 3PP during the control of our BCI using graph theory properties such as degree, betweenness centrality, and efficiency. Changes in these metrics are obtained for the case of the 1PP, as well as for the traditional third person perspective (3PP) in which the user can see the movement of the robot as feedback. As proof-of-concept, electroencephalography (EEG) signals were recorded from two subjects while they performed MI to control the movement of the robot. The graph theoretical analysis was applied to the binary directed networks obtained through the partial directed coherence (PDC). In our preliminary assessment we found that the efficiency in the α brain rhythm is greater in 1PP condition in comparison to the 3PP at the prefrontal cortex. Also, a stronger influence of signals measured at EEG channel C3 (primary motor cortex) to other regions was found in 1PP condition. Furthermore, our preliminary results seem to indicate that α and β brain rhythms have a high indegree at prefrontal cortex in 1PP condition, and this could be possibly related to the experience of sense of agency. Therefore, using the PDC combined with graph theory while controlling a telepresence robot in an immersive system may contribute to understand the organization and behavior of brain networks in these environments.


INTRODUCTION
A brain-computer interface (BCI) is a system that enables a real-time user-device communication pathway through different types of brain activity. In the beginning, BCIs were aimed for people with a disability of motor control or speech. Nowadays, even healthy people is using this technology for different applications. Some examples of BCI with immersive applications include gaming (Lalor et al., 2005), training (Vourvopoulos et al., 2016), rehabilitation (Calabrò et al., 2017), or psychological treatments (Jäncke et al., 2009). Furthermore, there have been some applications in which a user can teleoperate a robot, i.e., the user can have control of a robot that is not placed in the same location as him/her. Some examples of such applications are shown in Escolano et al. (2012), Leeb et al. (2015), Beraldo et al. (2018). In these cases, the subject perceives the environment real and in 3D as a extension of his/her sensorial functions. Such extension increases the feeling of presence of a remote scenario as well as a sense of agency when moving (Furht, 2008). This can be achieved by using either technologies based on head-mounted display (HMD) or multiple projections. Immersive virtual reality (VR) can also use HMD to project the virtual space just in front of the eyes, then the users focus on the display without distraction.
In the last years, there has been an increase in the number of research about sense of presence, embodiment, and sense of agency when combining BCI and immersive environments. Nevertheless, most studies of sense of presence have focused in the subjective experience analysis through questionnaires, like in Friedman et al. (2007) and Baka et al. (2015). Only a few studies involve the analysis of brain activity with the aim to explain cognitive processes related to the sense of agency in VR environments (Baumgartner et al., 2008). Therefore, the purpose of this study is to further contribute to this area with a quantitative analysis of the brain connectivity while controlling a BCI teleoperated robot in an immersive environment, in our case through graph theory metrics. For that purpose, we first describe the BCI which is controlled with motor imagery (MI) and two conditions: the first person perspective (1PP), i.e., the immersion experience, and the traditional third person perspective (3PP) of visual feedback. Next, we introduce the partial directed coherence (PDC) as the metric that allows us to assess the functional brain connectivity by calculating a connectivity matrix based on electroencephalography (EEG) measurements at various brain frequency bands: θ (4-7 Hz), α (8-13 Hz), β (14-29 Hz), and γ (30-50 Hz).
We already showed the usefulness of the PDC for the analysis of functional brain connectivity in Gaxiola-Tirado et al. (2018) and Alanís-Espinosa and Gutiérrez (2018). Yet, there is still work to do in how to interpret the interactions that the PDC reveals. The brain functional connectivity represents a complex system because of the transient nature of the interactions, such as synchronizations and desynchronizations between different brain regions. Therefore, it is complicated to compare the cortical connections within different brain rhythms. For this reason, here we propose the use of graph theory to describe the involvement of the different EEG channels (based on the PDC values) over the frontal, central, parietal, and occipital regions. Graph theory is a powerful tool for understanding the interactions and topology of various types of networks, and it has found a place in neurosciences for investigating brain aging (Vecchio et al., 2014), and different brain disorders like in mild cognitive impairment (Berlot et al., 2016), and brain function in spinal cord injured patients (de Vico Fallani et al., 2007). In BCI research, graph theory has been already proposed to analyze brain networks of different mental tasks (Huang et al., 2016;Stefano Filho et al., 2018). In this paper, we propose the use of a series of graph theory metrics in order to understand the differences in functional brain connectivity that arise depending on the immersion experience that our users have with our proposed BCI system.

MATERIALS AND METHODS
As proof-of-concept, we acquired the EEG recordings from two healthy right-handed participants (25 years-old female and 24 years-old male) that were asked to operate our proposed BCI system. Participants were recruited from the Center for Research and Advanced Studies at Monterrey and they were not familiar with BCI's. All experiments were conducted in an ethical and responsible manner, and with the approval of our institution's ethical committee. Informed consent was obtained from the participants after explanation of the study.
The experimental protocol consists of two stages. The first one is the conventional BCI training (2D, monitor-based), in which imaginary movements are performed according to the visual cues presented without feedback. The second stage corresponds to the BCI control task. Next, we will explain both of them.

Conventional BCI Training
The aim of this stage is to extract characteristic features from the EEG recordings to be used for the control of the BCI. First, the subject is required to either perform imaginary movements or to remain in resting state. Such behaviors are triggered by visual cues presented on a computer monitor while the subject is sitting on a chair. The cues are red arrows that indicate to continuously imagine an action based on the direction of the arrow: • or indicate the subject to open and close the left or right hand, respectively; • indicates to open and close both hands; • indicates to move up and down both feet; • If there is not a red arrow in the screen, then the subject should remain at rest.
Hence, the subjects are tested for four different control conditions or classes defined by the arrows. Each trial starts with a fixation cross in which the subject should be at rest. After 3 s of fixation, the visual cue is presented for 4 s. After that, a black screen is shown indicating the subject to rest. Each trial is separated by a random duration interval of 2.1-2.5 s between them to avoid adaptation (Friedman et al., 2007). The training sequence is shown in Figure 1. The subject has five minutes of rest between each one of the five runs. Each run corresponds to the acquisition of EEG measurements of forty trials (ten randomized presentations of the cues of each class). These trials are used as calibration recordings to build a classifier for discriminating two different mental tasks (in our case, each of the tested classes compared against the resting state). The task that shows highest separability vs. resting state is then used as the user-specific input for the BCI control task system in the next stage of the experiments.
The EEG signals acquired under the previously explained protocol are filtered and analyzed with BCI2000 offline tools to get a personalized feature for BCI control task. This is done by the calculation of the coefficient of determination for all electrodes within 0-70 Hz for two cognitive states. The coefficient of determination, or r 2 value, is a statistical measure computed over a pair of sample distributions, and it provides a measure of how strongly the means of the two distributions differ in relation to variance. In a BCI context, r 2 values are computed over signals that have been measured under two different task conditions, and they represent the fraction of the total signal variance that is determined by the task condition. In this way, the channels and frequencies with higher r 2 values are selected as features to train the classifier. In Figure 2, an example of the results for the power spectra and r 2 values of the channels used to train the classifier are shown. Note that the power in C3 at 11 Hz corresponding to MI of right hand is lower in magnitude than in the resting condition, and this is expected because of the desynchronization of µ rhythm at the motor cortex. More details on the process of r 2 calculation can be found at http://www.bci2000.org/mediawiki/ index.php/User_Tutorial:Mu_Rhythm_BCI_Tutorial.

BCI Control Task
For this stage, we used the NAO humanoid robot as the device the BCI system is to control. The main goal of the subject is to control the movement of the robot with the personalized classifier obtained in the training stage. The sequence displayed on the  screen to subjects during the experiment is shown in Figure 3. At the beginning of the experiment, the subject remains in resting state for 15 s. Then, the following sequence is displayed on the screen: it starts with a green cross, followed by a black screen or a red arrow, and finally the letters FB (for feedback) appear. The first and second screen appear for 4 s, and the FB screen appears for 10 s with a random time between each trial from 100 to 500 ms with a black screen. During the FB screen, the subject receives feedback from the robot moving or not according to the cue. When a MI task is detected, the robot closes and opens its hand, otherwise, the robot does not move when resting state is detected. The two types of feedback implemented are the following: 3PP: The subject is comfortably sitting in front of the robot and a monitor where the sequence of stimuli is shown, as it can be seen in Figure 4A. In this case, the subject sees the robot moving in third-person perspective. For all trials, the robot is placed on the left side, and the monitor on the right side.
1PP: The subject is comfortably sitting using a passive HMD in which it is shown whatever the robot sees, as well as the stimuli sequence (see Figure 4B). Meanwhile, the robot is placed outside the room as shown in Figure 4C. The implementation of this feedback scheme is further detailed in section 2.3.
The subject is instructed to avoid blinking during the MI task to minimize noise in EEG data. Additionally, the subject has 5 min of rest between runs. A total of 200 trials of 18 s of duration were recorded in two different sessions per perspective. The accuracy is recorded for each cue, which corresponds to the number of times that the subject correctly controlled the movement of the robot, as well as the number of times the robot halted during the subject's resting state.

Immersive Telepresence System
The aim of this system is to allow the BCI user to see the world through the perspective of the NAO robot, i.e., as if the robot's movements were the user's own. In order to be immersive, some requirements have to be considered in the design of the telepresence: 1. BCI user must control the NAO robot placed in a remote site using the MI paradigm; 2. control signals have to be sent wirelessly; 3. control signals and video feedback must have the minimum lag for an immersive experience; 4. BCI user has a stereoscopic video image as feedback displayed into a HMD.
The implementation of the immersive telepresence system can be divided into two parts. The first contains the components of the BCI system considering the stereoscopic video feedback from the remote environment. The second part comprises the BCI software implementation, the system for the communication channel, and the system implemented at the remote environment. The general setup can be seen in Figure 5.
The hardware of the immersive telepresence system can be divided in three different modules: • Module I corresponds to the location of both the BCI system and its user. It consists of a desktop computer (Intel core i3, 8 GB RAM) and the signal acquisition system MOBITA, which is a wireless EEG system with 32-channel. These two pieces are connected via 2.4 GHz wireless network. The virtual reality headset consisted of a passive HMD, coupled with a Samsung Galaxy S6 cellphone. The Galaxy S6 has a large screen (1,440 × 2,560), an Exynos processor 7420 2.1 GHz, 3 GB of RAM and GPU Mali-T760MP8, and it is powerful enough to have an unnoticeable delay of the video streaming and VR web application which allows the immersion feeling. • In module II, at the NAO's robot site, a helmet with an Arducam of 5 MP and 1080p video resolution is placed on top of the robot's head in order to provide the required video feedback that the native robot's webcam cannot. The Arducam is connected directly to Raspberry Pi's native CSI camera port to provide better performance than a webcam in terms of the frame rate and resolution. • Module III corresponds to the communication channel, which is implemented through an ASUS AC1200 router of double band and links to the main server running in the RPi 3.
The software architecture is shown in Figure 6. In our case, the software used at the BCI user site is running in a desktop computer with Microsoft Windows 7 of 64 bits as operating system which runs Python 2.7 and OpenViBE (Renard et al., 2010). Within OpenViBE, the Python box is used to send the control signals to the RPi 3 and they are echoed to control robot using TCP/IP. Meanwhile, the stimulus provided by a LUA stimulator box is sent using an HTTP server written in Python (see www.lua.org for more details about LUA scripting language). The stimulus is received with a HTTP client at the RPi 3 and echoed to the main server. At the teleoperator site, the RPi Cam-Web-Interface was implemented to manage the streaming video from the mini camera. Such interface is freely available at https:// github.com/silvanmelchior/RPi_Cam_Web_Interface. In order to obtain a stereoscopic video image, we decided to use only the video from one camera as in the Pi Viewer demo available at https://github.com/patcat/PiViewer, in which the video from the RPi camera is streamed into VR using JavaScript. In our case, it was implemented with some modifications to enable the appearance of the stimulus images from OpenViBE without a noticeable delay into a web app.
The main server running in the RPi 3 was written with Node.js (https://nodejs.org). This server allows the integration of the   Frontiers in Psychology | www.frontiersin.org camera and stimulus into a web app. It hosts the web page where it is shown a stereoscopic image of the remote environment streamed by the mini camera, as well as the different stimuli synchronized from the BCI user site, as it is shown in Figure 7.

EEG Data Acquisition and Preprocessing
Subjects were comfortably seated on a chair inside a noise free and normally lighted room. A 32-channel EEG system (Mobita TMSi) was used to record the brain electrical potentials by means of an electrode cap with sensors placed according to the 10-20 international system and with reference to AFz. Impedance of all electrodes were kept below 5 k . The acquisition was performed at a sampling rate of 1,000 Hz, the signals were bandpass filtered with a zero-phase fourth-order Butterworth between 1 and 100 Hz band and a notch filter to remove artifact caused by electrical power lines in 60 Hz. The blinking artifacts were removed using independent component analysis (ICA) and finally a baseline correction was performed.

Feature Extraction and Classification
The EEG signals obtained from the subjects after the conventional BCI training scheme described in section 2.1 were analyzed with the general-purpose software system BCI2000. This software allows for the calculation of personalized features for BCI control, and this is done by the calculation of the r 2 values for all electrodes within 0-70 Hz for two cognitive states. In this way, the channels and frequencies with higher r 2 values are selected as features to train the classifier. More details on this process can be found at Mu Rhythm BCI Tutorial.
Once the features have been extracted, a linear discriminant analysis (LDA) classifier (Kantardzic, 2002) is used to discriminate between the control command (detection of features that characterize the MI task) or the rest state.

Partial Directed Coherence
The partial directed coherence (PDC) is a method based on the Granger causality to measure the coupling or connectivity between different channels in the frequency domain. The PDC also identifies the direction of information flow and the strength (Baccala and Sameshima, 2001). In our case, the FIGURE 7 | Stereoscopic image and BCI stimulus. The blue arrow is the stimulus to perform the motor imagery movement. frequency domains selected correspond to the known EEG rhythms θ , α, β, and γ . The PDC was calculated only between the following EEG channels: FP1, FP2, F3, F4, Fz, C4, C3, Cz, C4, C3, Pz, T5, T6, O1, Oz, and O2. We chose those channels based on previous work related to brain processing in the human visuomotor system (Binkofski and Buxbaum, 2013), and during motor imagery events (Ghosh et al., 2015).
In order to compute the PDC, the data needs to be fitted to a multivariate autoregressive (MVAR) model. For the case of a set S = {x m , 1 ≤ m ≤ M} of M EEG signals (in our case those from the previously selected M = 16 channels), the MVAR model of order ρ of x(n) = [x 1 (n), x 2 (n), . . . , x i (n), . . . , x j (n), . . . , x M (n)] T , for n = 1, 2, . . . , T time samples, is given by where A 1 , A 2 , . . . , A ρ are the M × M coefficient matrices containing the coefficients a ij (l) accounting for the linear interaction effect of x j (n − l) onto x i (n), and v(n) = [v 1 (n), v 2 (n), . . . , v M (n)] T is the noise vector (uncorrelated error process). Under those conditions, a measure of the direct causal relations (directional connectivity) of x j to x i for a frequency f is given by the PDC as where A ij (f ) and a j (f ) are, respectively, the i, j element and the j-th column of In this work, the connectivity analysis is done during MI epoch from the first to the third second after the cue is presented in

P2, 1PP
FIGURE 9 | Significant connectivity (left) and its corresponding degree distributions (right) for the 3PP (A,C) and 1PP (B,D) conditions in α band for P1 and P2, respectively. Indegree distribution is represented by dark-color bars and the outdegree is represented by light-color bars.
order to calculate the PDC values for each frequency between 4 and 50 Hz. All values are arranged in matrices of size 16 × 16 (to account for all π ij values of our EEG channels) for each of the frequencies. Then, the maximum PDC value is calculated between those PDC matrices corresponding to the frequencies of each brain rhythms. At the end, we are left with four PDC matrices per trial (representing θ , α, β, and γ ), which are then used for the graph analysis.

Graph Analysis
Graph theory has been used to describe large scale networks in different research fields. In neuroscience, graph theory is employed as a network analysis to identify the simultaneous activity of different brain regions stimulated by a mental state. In our case, a brain network is obtained based on the PDC matrix generated with the PDC values of each pair of channels {i, j}. From the perspective of graph theory, the set of nodes correspond to the EEG channels, while the set of edges represent an anatomical link between those nodes or a functional dependence (Sporns, 2010). These pairwise connections are accounted for in a connection matrix.
In order to compute the topological features of interest, the connectivity matrix is converted into a directed, unweighted graph (also referred to as digraph). Such unweighted matrix is then binarized by choosing a threshold τ that represents the number of most powerful connections (connection density). In our study, τ represents the ratio between the effective connections and all the possible ones in the digraph. A range of τ = [0.3, 0.9] was explored for each frequency band with increments of 0.05. Then, we identified the minimum p-value in which there was a significant difference (p ≤ 0.05) of the global efficiency metric when comparing it between the 1PP and 3PP cases. Figure 8 shows an example of such search when p-values are obtained from a Wilcoxon rank-sum test for the case of τ between 1PP and 3PP in α band. Table 1 summarizes the results for τ in each frequency band, which then we used to obtain the digraphs for the following graph metrics.
Once it was binarized, the digraph was characterized according to the degree, distribution degree, node betweenness centrality, and local efficiency. All those network metrics were computed using the brain connectivity toolbox for Matlab (Rubinov and Sporns, 2010). Next, a more detailed description of each of them is presented.

Degree Distribution
The degree k of a node i measures the number of connections, so it indicates how many nodes are connected to node i. The degree distribution P(k) represents the probability to find a node i with certain degree k. In our case, the degree distribution can be divided in two: the indegree (ID) denoted as k in , and the outdegree (OD) denoted as k out . Those represent the total number of connections incoming to a node and outgoing from a node, respectively. A large value for k in means that the node is influenced by a large number of different channels. A large value of k out means a node has a large number of potential targets.

Betweenness Centrality
In a graph structure, we can identify important nodes that often interact with many others as a way to facilitate functional integration (Rubinov and Sporns, 2010). Within a graph, we can also find central nodes that participate in many short paths, and these act as important controllers of information flow. Betweenness centrality is a sensitive measure defined as the fraction of all shortest paths in the network that pass through a given node. This measure is used to identify important functional connections.

Efficiency
This metric was introduced by Latora and Marchiori (2001) to measure how efficiently the nodes communicate between them. The efficiency e ij of the communication between the nodes i and j is inversely proportional to the shortest path length d ij . If a path does not exist between the nodes i and j, then d ij = ∞ and e ij = 0. The global efficiency is defined as  where N corresponds to the number of nodes composing the graph.

Statistic Assessment of Connectivity
Based on eConnectome toolbox (He et al., 2011), a nonparametric method based on surrogate data is used to assess the statistical significance of our calculated values of PDC. The Fourier transform is applied to each trial of the original data to randomly shuffle the phases without changing the magnitude as follows. Next, inverse Fourier transform is applied to each trial with permutated phase to generate surrogate data in time, followed by the PDC estimation. The shuffling and PDC calculation are repeated 1,000 times, resulting in a distribution of PDC values under the null hypothesis that no connectivity exists and with a significance level of 0.05 (Yasumasa Takahashi et al., 2007).

RESULTS
The analysis within the two types of experiments was done with the MI of the right hand because it was the one that presented the highest values of r 2 for Participant 1 (P1) and Participant 2 (P2). Therefore, the movement of the NAO robot was chosen to be also the closing and opening of its right hand to promote the sense of agency.
Next, the degree distributions for ID and OD were calculated for each frequency band. As an example, the significant connectivities and their corresponding degree distributions in α band for both the 1PP and 3PP conditions are shown in Figure 9 for P1 and P2, respectively. Histograms are normalized to the number of possible nodes in the network (15 nodes). We note the consistency in the behavior of the distribution in both participants. An interesting result, is that the OD distribution shows different trends within each condition. The right-skewed tails of the OD distributions indicate that there are no nodes with less than 7 outgoing links. Contrary to the ID distribution, in which there are no nodes with more than 7 incoming connections.
Similarly, Figure 10 shows the comparison between 1PP and 3PP of the ID in β for all the channels of interest of each participant. We can identify in which node a specific metric is greater or lower for both participants. For example, we can identify a greater ID in 1PP than in 3PP at FP1, which would mean that FP1 is dependent by more channels in 1PP than in 3PP. Likewise, Oz has a greater ID in 3PP than in 1PP. This means that the dependency of channel Oz from other channels is lower in 1PP than in 3PP.
In Figures 11-14, we show a summary of all the graph metrics obtained from our participants in the different brain rhythms. These results are presented over the array of sensors for a better  view of the brain regions they correspond to. Furthermore, only the metrics that increased in both participants are shown for both 1PP and 3PP conditions: • Figure 11 shows the magnitude of the indegree. We can identify the stronger value of the ID in β, θ , and γ for 1PP condition at right temporal, parietal, and occipital brain regions. • Figure 12, shows the magnitude of the outdegree. We can identify in β that at C3 and the midline it is greater in 1PP than in 3PP. This could mean that the influence of these channels is stronger than the others. Additionally, there is a high OD at T6 in the 3PP condition in all frequency bands, which could mean that this channel is an important hub in 3PP condition. • Figure 13 shows the magnitude of the node betweenness centrality. The greatest betweenness centrality is at C3 in α for the 1PP condition. In θ and γ the betweenness centrality nodes are at the temporal and occipital regions in 1PP condition. Meanwhile, for the 3PP condition are at the frontal region. We didn't find any significant difference in α. • Figure 14 shows the magnitude of the local efficiency. We can identify high efficiency in β for both conditions, but for 1PP is greater at P3, while from 3PP is greater at C3. We did not find any significant differences in the other frequency bands. As an important result, we identified a greater efficiency in 1PP than in 3PP at the prefrontal area (FP1 and FP2). This could be due to the sense of agency, similarly to Baumgartner et al. (2008) where they found the prefrontal areas strongly involved in the modulation of experience of presence in adults.

DISCUSSION
This proof-of-concept study, whose main purpose is to showcase an experimental platform with great potential to facilitate the user's immersion experience, also allowed to exemplify the use of graph theory to investigate relevant features in brain networks through EEG data. Such data was acquired from an experimental BCI system that is well suited to describe the differences of sense of agency in two environments with different levels of immersion. We calculated functional connectivity based on PDC because it is a metric that allows for an analysis in the frequency domain, then giving us more information about the cognitive processes in a specific band. This brain network analysis could help us understand the causal relationships and give us an idea of the information flow and brain organization during the control of an immersive BCI and sense of agency.
Our results show the frontal and parietal brain regions as target areas during 3PP, mainly in β. This is consistent with Rathee et al. (2016) who reported that, using the partial Granger causality during the right hand MI, the higher FIGURE 13 | Mean node betweenness centrality for both participants shown in color only at the sensors with an incremental trend in comparison to the other condition. For each frequency band, the heads to the left (A,C,E) show the cases when the node betweenness centrality in 1PP increased in comparison to 3PP, while the heads to the right (B,D,F) show those for which 3PP increased in comparison to 1PP. No increase in node betweenness was found in α. FIGURE 14 | Mean efficiency for both participants shown in color only at the sensors with an incremental trend in comparison to the other condition. The first two heads (A,B) show the cases when the efficiency in 1PP increased in comparison to 3PP, while the third head (C) shows those for which 3PP increased in comparison to 1PP. No increase in node efficiency was found in α for 3PP, nor in θ and γ in any modality. amount of outgoing information was from the central region (C3, Cz, C4) and the targets were Fz and Pz. Furthermore, the nodes with higher ID in 3PP in α and β are in accordance with the results in Athanasiou et al. (2018), where the nodes with higher in-degree are associated with the supplementary motor areas bilaterally during hand motor imagery, which corresponds to F3, Fz, and F4. Moreover, α and β presents a high ID at FP1 in 1PP condition, which might possibly be related to a better coordinating role of these rhythms in the planning of the sensorimotor process. Also, Baumgartner et al. (2008) found that the prefrontal areas are strongly involved in modulating the experience of presence. In Ghosh et al. (2015), they found as characteristic the out strength of C3 while performing motor imagery of the right hand. In our case, we found that the OD at C3 in 1PP is higher than in 3PP. Then, the stronger influence from C3 to other regions might be possibly a reflection of a better engagement.

Cz
In this study we found a greater OD at the SMA (C3, Cz) in 1PP than in 3PP in β and γ , which could be related to the regulative role of the SMAs during the planning of a movement, as identified in Athanasiou et al. (2012). Furthermore, when performing MI of right hand in 1PP, there is a greater OD and node betweenness centrality in β than in 3PP, which can be an indication of the existence of a central node for information exchange. This is in line with results from Li et al. (2019).
The node betweenness centrality shows us an important functional connection at C3 in β, which could mean that C3 plays the role of regulator in the flow of information during the MI in 1PP condition. As seen in α, the efficiency is greater in 1PP condition than in 3PP at the prefrontal cortex, which could be considered as a structural candidate for modulating inter-individual differences in the experience of presence (Baumgartner et al., 2008).
The present study is limited by its small sample size in human subjects, but our primary goal was the development and demonstration of our experimental BCI system and the analysis approach of the data we collected from it. Future work would include a more exhaustive experimentation in different subjects, as well as increasing the sense of agency by incorporating other activities and other interfaces, such as the one described by Gaxiola-Tirado et al. (2019). Additionally, and given that in this work we only analyzed the data during the MI epochs, we would like to extend the analysis to the feedback epochs as well so to discard the possible involvement of mirror system neurons (Pineda, 2005) and properly assess the experience of sense of agency.

CONCLUDING REMARKS
In this paper, we showed the applicability of an experimental BCI system for the interaction with a robot through different immersion experiences. The BCI system we proposed can be seen as a tool for the analysis of functional brain connectivity associated to the control of a BCI in an immersive environment. For that purpose, we expended our previous work on connectivity measures based on the PDC by using metrics of graph theory. As a proof-of-concept, we analyzed the data from two participants, and the results seem to be in line with previous results from the literature.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Internal Committee of Ethics, Cinvestav Monterrey. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
MA-E and DG conceived and designed the study and proofread and edited the paper. MA-E collected the data, performed the signal processing and statistical analysis, contributed to the interpretation of the findings, and drafted the paper. All authors read and approved the final version of the paper.

FUNDING
The work of MA-E was sponsored by the Mexican Council of Science and Technology (CONACyT) through the graduate school scholarship no. 222676.