Neural Plasticity
Volume 2018 (2018), Article ID 9354207, 20 pages
https://doi.org/10.1155/2018/9354207

Clinical Study

Functional Brain Connectivity during Multiple Motor Imagery Tasks in Spinal Cord Injury

Alkinoos Athanasiou,1,2 Nikos Terzopoulos,¹ Niki Pandria,¹ Ioannis Xygonakis,¹ Nicolas Foroglou,² Konstantinos Polyzoidis,² and Panagiotis D. Bamidis¹

¹Biomedical Electronics Robotics & Devices (BERD) Group, Lab of Medical Physics, School of Medicine, Faculty of Health Sciences, Aristotle University of Thessaloniki (AUTH), 54124 Thessaloniki, Greece
²1st Department of Neurosurgery, “AHEPA” University General Hospital, Aristotle University of Thessaloniki (AUTH), 54636 Thessaloniki, Greece

Correspondence should be addressed to Alkinoos Athanasiou; athalkinoos@auth.gr

Received 29 December 2017; Revised 6 March 2018; Accepted 21 March 2018; Published 2 May 2018

Academic Editor: Thierry Pozzo

Copyright © 2018 Alkinoos Athanasiou et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Reciprocal communication of the central and peripheral nervous systems is compromised during spinal cord injury due to neurotrauma of ascending and descending pathways. Changes in brain organization after spinal cord injury have been associated with differences in prognosis. Changes in functional connectivity may also serve as injury biomarkers. Most studies on functional connectivity have focused on chronic complete injury or resting-state condition. In our study, ten right-handed patients with incomplete spinal cord injury and ten age- and gender-matched healthy controls performed multiple visual motor imagery tasks of upper extremities and walking under high-resolution electroencephalography recording. Directed transfer function was used to study connectivity at the cortical source space between sensorimotor nodes. Chronic disruption of reciprocal communication in incomplete injury could result in permanent significant decrease of connectivity in a subset of the sensorimotor network, regardless of positive or negative neurological outcome. Cingulate motor areas consistently contributed the larger outflow (right) and received the higher inflow (left) among all nodes, across all motor imagery categories, in both groups. Injured subjects had higher outflow from left cingulate than healthy subjects and higher inflow in right cingulate than healthy subjects. Alpha networks were less dense, showing less integration and more segregation than beta networks. Spinal cord injury patients showed signs of increased local processing as adaptive mechanism. This trial is registered with NCT02443558.

1. Introduction

Reciprocal communication of the central and peripheral nervous systems is compromised during spinal cord injury (SCI), a condition that often causes permanent disability due to massive neurotrauma of ascending and descending pathways [1, 2]. While changes in brain activity and brain organization may seem trivial, when compared to the underlying injury of the pathways, they have nevertheless been consistently associated with SCI [3–5]. Such changes have also been observed at the early stages after the injury and have been associated with differences regarding the prognosis of SCI patients’ recovery [4–6]. Demonstrated structural changes of the brain include atrophy of afferent neural pathways, microstructural changes of efferent axons, and disorder of white matter integrity in multiple nodes of the sensorimotor cortex that involve the primary motor and somatosensory areas [7] and also diffuse neuronal degeneration [8].

Functional connectivity (FC) after SCI has been studied by means of electroencephalography (EEG) [9–15] and functional magnetic resonance imaging (fMRI) [6, 16–21]. Poor recovery after SCI has been associated with decreased FC strengths between midline sensorimotor network nodes during resting state, while the opposite pattern has been associated with good recovery [6]. Supplementary and cingulate motor areas have been shown to play important roles during the sensorimotor neurophysiological process [9, 11], while unique interactions and temporal dynamics have been identified in the functional networks of SCI patients [12, 14]. Connectivity changes have been hypothesized to be able to serve as injury biomarkers [22] while novel methods have been developed and have been proposed in order to study the brain’s connectome following SCI in hopes of providing more reliable evidence of these changes [23].

Most studies on functional connectivity after SCI have focused on patients with chronic complete injury, including the majority of EEG studies. Only a couple of studies employing fMRI as a modality to detect brain activity have assessed patients with incomplete injury but those have only studied connectivity during resting-state condition so far [6, 16–18]. Moreover, the pioneer EEG studies on functional cortical connectivity of SCI patients during motor imagery have employed a robust but rather limited study design [9]. In this design, the motor imagery task involved an attempt to move the paralyzed right foot and was performed simultaneously with one motor execution task (lip protrusion), while the functional networks were subsequently analyzed using a variety of tools. So far, to the authors’ best knowledge, no study has been performed employing multiple motor imagery tasks, especially of the upper extremities, aiming to analyze differences in the formed networks. Moreover, incomplete injury at the chronic phases remains understudied compared to chronic complete injury with regard to functional connectivity networks. Despite the clinical and social impact of SCI, so far published studies have been unable to form a complete model regarding the effect of the injury on brain networks, although effort has been made into modeling-specific aspects, like resting-state connectivity and chronic injury [3, 22]. It can be hypothesized—and there are also indications—that even incomplete spinal cord injury may show measurable effects on the functional sensorimotor network [24] that could be important for modeling the condition in relation to prognosis [3, 6, 16].

Motor imagery, apart from its importance to the study of brain activity after neurotrauma, has shown great potential in motor skill learning and in rehabilitation of upper and lower limb paralysis [25, 26]. It has been established that motor imagery produces patterns of brain activation and brain connectivity similar to those of motor execution [27, 28], while the visual motor imagery class also activates a distinct task-dependent neural system [29, 30]. Motor imagery has been used as a modality to induce plasticity and recovery in a range of conditions [31], including complete cervical spinal cord injury [32] and stroke [33]. Moreover, motor imagery has been also used as a control modality for brain-computer interface implementations of exoskeletons for complete [34] and incomplete spinal cord injury [35]. Functional recovery has been induced even in the case of complete injuries using such an approach [34]. Such results demonstrate the importance of motor imagery functional networks studies to accurately model the plastic changes that occur after SCI.

In our previous work, we have presented our study with a cohort of SCI patients and healthy control subjects that exercised motor imagery to achieve control of anthropomorphic robotic arms in various movement tasks. We have accounted for development [36], pilot experiments, and brain network analysis [37, 38], and we have presented a detailed user perception and performance assessment study, based on neurological and psychometric evaluation [39]. In the current paper, we present an elaborate analysis of the functional connectivity networks formed on the sensorimotor cortex during visual motor imagery of multiple motor tasks performed by subjects with SCI and healthy controls. In Materials and Methods, we briefly present the experimental setup and detail our signal processing computational workflow, network analysis, and statistical comparisons. In Results, we detail important findings with regard to the effect of injury, motor imagery category, brainwave rhythm, and timing of imagery. In Discussion, we attempt to interpret our results in the context of already published studies in the field, and we note the limitations of this approach.

2. Materials and Methods

2.1. Experimental Setup

2.1.1. Recruitment and Subject Assessment

The experimental setup has been previously described in detail, including subject assessment [39] and procedures [36, 37], so we will hereby provide only a brief overview. Our experimental protocol was approved by the institutional ethical committee [40], and all subjects signed an informed consent form. We recruited 8 male and 2 female patients with SCI (age: mean 46.0, SD 17.64, range 28–74 years) and ten gender- and age-matched healthy controls. All participants were right-handed and reported no prior experience in mental imagery (Table 1).

Table 1: Demographic data and reported imagery capacity of subject groups (SCI and healthy).

(a)

SCI group
Age
Gender
Cause
AIS
NLI
VVIQ

CSI-02-001
28
f
MVA
C
C4
54

CSI-02-002
52
m
MVA
D
C4
69

CSI-02-003
42
m
MVA
D
C8
68

CSI-02-004
70
m
Fall
D
C5
76

CSI-02-005
60
m
Fall
E
C6
70

CSI-02-006
28
m
MVA
D
C5
56

CSI-02-007
30
m
MVA
E
C5
67

CSI-03-001
47
m
Fall
A
T7
72

CSI-03-002
29
f
MVA
B
T4
60

CSI-03-003
74
m
Other
B
T4
65

Mean (SD)
46.00 (17.64)

65.70 (7.04)

(b)

Healthy group
Age
Gender
VVIQ

CSI-04-001
27
f
77

CSI-04-007
51
m
75

CSI-04-003
43
m
56

CSI-04-006
71
m
70

CSI-04-009
63
m
70

CSI-04-004
28
m
46

CSI-04-005
31
m
58

CSI-04-008
47
m
80

CSI-04-002
27
f
75

CSI-04-010
74
m
63

Mean (SD)
46.20 (18.27)

67.00 (10.09)

For both groups, we collected demographics and medical history; also, a specialist physician performed neurological examination using the International Standards for Neurological Classification of Spinal Cord Injury, to document classification in American Spinal Injury Association Impairment Scale (AIS) and the Neurological Level of Injury (NLI). Subject assessment also included subjective reporting of imagery capacity, using Vividness of Visual Imagery Questionnaire (VVIQ) [41] with eyes open (Table 1). Within the SCI group, 60% of the patients were grouped into positive outcome based on the neurological assessment (4 AIS D, 2 AIS E), and 40% of the patients were grouped into negative outcome (1 AIS A, 2 AIS B, and 1 AIS C).

2.1.2. Experimental Procedure

The experiment took place inside a magnetic shielded room for EEG recording, specially designed for presentation capability and audiovisual monitoring of the participants. The subjects sat at a 1 m distance across a 21 computer monitor. They wore an active electrode cap (Brain Products, Germany) and were connected to a 128-channel EEG (Nihon-Kohden, Japan) according to the high-resolution EEG 10–5 international electrode system [42]. Recordings were taken at a sampling rate of 1000 Hz and an impedance threshold of 10 kohm [36]. Initially, we recorded resting-state activity, 1.5 min with open eyes and 1.5 min with closed eyes.

In the main experimental part, the subjects watched videos of 32 different arm motor tasks, a walking task video, and an oddball video. The videos were presented in random order. All arm motor task videos were presented from the perspective of the participant watching his or her own arms and were gender-matched. The walking task video presented a pair of gender-matched legs walking, while seeing them from the perspective of watching one’s own legs. The oddball video showed a wildlife documentary. The videos lasted 5 sec, each followed by a black screen with duration of 4 sec. A trigger channel was recorded at the onset of each visual cue (the start of each video) through an optic diode. The presentation was separated in three parts of about 17 min each, with 5 min of rest between them (Figure 1). At conclusion, each task had been presented 9 times in total. The subjects were asked to perform visual motor imagery (VMI), while watching a motor task, without actually moving their limbs (regardless of neurological status or group) and were instructed to rest while watching the black screen. They already knew that the videos would be presented at random but not of the presence of an oddball video (video showing wildlife). Moreover, the subjects’ arms, torso, and legs were covered with a black curtain during the whole procedure in order to facilitate mental registration of the presented arms and legs into their perceived body schema [43].

Figure 1: Flow diagram of the experimental procedure of one part of the visual motor imagery presentation. Each presented video lasted 5 seconds and was followed by 4 seconds of black resting screen. The videos were presented 9 times each, in a random order. The presentation was divided into three parts, lasting approximately 17 minutes each, with an intermission between them.

The 32 upper extremity motor tasks consisted of 8 independent movements (degrees of freedom) 2 directions of movement 2 extremities, comprising the full range of motion of the human arms and were classified into categories for further analysis [39]. In short, the 8 categories of motor tasks were “Hands,” “Left,” “Right,” “Proximal,” “Distal,” “Rotational,” “Linear,” and “Walking.”

The “Hands” category included all 32 tasks of both upper extremities. The “Left” and “Right” categories included 16 motor tasks each of the respective upper extremity (left or right). The “Proximal” category included 16 motor tasks of the shoulder and elbow joints of both extremities, while the “Distal” category included the remaining 16 motor tasks of wrist joints and fingers. Further, the “Rotational” category included those 8 motor tasks that result in rotational motion, and the “Linear” category included those 24 motor tasks of both extremities that resulted in linear motion (Table 2). Finally, the “Walking” category was also defined as a separate category of motor imagery, consisting only of the walking motor imagery task, for a total of 8 categories. Table 2 lists in summary all presented motor tasks of one upper extremity (16 tasks), for each showing the degree of freedom, direction of movement, classification by proximity, and resulted motion.

Table 2: Presented motor tasks for one upper extremity (left or right): 16 motor tasks were presented (8 independent movements (degrees of freedom) 2 directions of movement) and were then classified by proximity (proximal or distal tasks) and resulting motion (linear or rotational). For both upper extremities, the subjects watched and performed visual imagery of 32 motor tasks in total.

Independent movement
Direction
Proximal/distal
Linear/rotational

Shoulder
Arm down
Proximal
Linear

Shoulder
Arm up
Proximal
Linear

Shoulder
Arm left
Proximal
Linear

Shoulder
Arm right
Proximal
Linear

Elbow
Forearm down
Proximal
Linear

Elbow
Forearm up
Proximal
Linear

Forearm
External rotation
Proximal
Rotational

Forearm
Internal rotation
Proximal
Rotational

Wrist
Hand down
Distal
Linear

Wrist
Hand up
Distal
Linear

Wrist
External rotation
Distal
Rotational

Wrist
Internal rotation
Distal
Rotational

Thumb
Open
Distal
Linear

Thumb
Close
Distal
Linear

Fingers
Open
Distal
Linear

Fingers
Close
Distal
Linear

2.2. Signal Analysis

2.2.1. Signal Preprocessing

Signal analysis was performed on a subset of the 10–5 international electrode system that is overlying the cortical sensorimotor areas [44] that were later defined as regions of interest (ROIs) for this study (Figure 2). This subset included 64 electrodes: AFF5h, AFF3h, AFF1h, AFz, AFF2h, AFF4h, AFF6h, F5, F3, F1, Fz, F2, F4, F6, FFT7h, FFC3h, FFC1h, FFC2h, FFC4h, FFT8h, FT7, FC5, FC3, FC1, FC2, FC4, FC6, FT8, FTT7h, FCC5h, FCC3h, FCC1h, FCC4h, FCC6h, FTT8h, C5, C3, C1, Cz, C2, C4, C6, CCP3h, CCP1h, CCP2h, CCP4h, CP5, CP3, CP1, CPz, CP2, CP4, CP6, CPP3h, CPP1h, CPP2h, CPP4h, P3, P1, Pz, P2, P4, PPO1h, and PPO2h. As scalp electrodes capture mixed activity from unknown cortical and subcortical brain sources, recording brain activity related to motor tasks only from the sensors overlying the sensorimotor area presents some risk for loss of information but also presents certain advantages. This approach has been used in EEG source imaging studies regarding motor tasks with good results [45–48], as the signal of interest is less attenuated and signal to noise ratio is higher than in distant sensors, while contaminated channels closer to muscular artifact generators are excluded.

Figure 2: Regions of interest (ROIs) of the sensorimotor cortex and the overlying subset of electrodes that was used for signal analysis in our study. 1: presupplementary motor area (pSMA); 2: supplementary motor area (SMA); 3: dorsal premotor area (PMd); 4: ventral premotor area (PMv); 5: cingulate motor area (CMA); 6: primary foot motor area (M1F); 7: primary hand motor area (M1H); 8: primary lip motor area (M1L); 9: primary foot somatosensory area (S1F); 10: primary hand somatosensory area (S1H); 11: secondary somatosensory area (S2); 12: somatosensory association area (SAC).

All signal preprocessing was performed using a custom script on the FieldTrip toolbox for MATLAB [49]. Raw data from those selected channels was band-pass filtered at 0.5–30 Hz using a zero-phase FIR filter and subsequently downsampled at 100 Hz. We visually examined continuous EEG signal time series of each subject to detect bad electrodes that showed large drifts from their mean value and then removed these electrodes. Epochs were then initially extracted from −2000 msec to +5000 msec centered on the trigger (visual cue). Subsequently, independent component analysis was performed on the concatenated continuous data (of each session) using the second-order blind identification method [50]. Independent components corresponding to eye blinks and muscle artifacts were identified and removed from the epoched data. Bad electrodes were then interpolated using spherical splines interpolation [51]. Finally, the epoched data were split into two time intervals (Figure 3), which will be referred to as “early” (early onset imagery from −1000 msec to +2000 msec around the trigger) and “late” (late continuous imagery from +2000 msec to +5000 msec after the trigger). While shorter time windows have also been used in similar analyses [52], differences in the behavior of alpha and beta rhythms between the window around the imagery onset and later windows have been identified with regard to networks [53], relative power [54], and event-related desynchronization [55]. The data from one subject (from the healthy group) was exempted from further analysis, as this preprocessing methodology did not result in sufficiently clean epoched data.

Figure 3: Activation time series of all regions of interest (ROIs) of the sensorimotor cortex during a random epoch and definition of time intervals around the trigger (onset of the video presenting the motor imagery task).

2.2.2. Current Cortical Density

The solution of the forward problem, the lead field matrix that best describes the conduction from the current cortical density (CCD) source model (Table 3) to scalp potentials, is based on the following equation:

where refers to the M simultaneous electrode voltage recordings, refers to the N current dipoles in the current cortical density model, is the noise vector, and is the abovementioned lead field matrix [56]. We used the solution applied in eConnectome toolbox for MATLAB [57, 58] of the forward problem which is a high-resolution lead field matrix relating 2054 scalp triangles to 7850 cortical dipoles. The lead field matrix is derived using a three-layer block element modifier model based on the Colin27 Montreal Neurological Institute brain [59]. The dipoles were constrained to the gray matter with orientations perpendicular to the local cortical surface, under the assumption that the primary source of measured EEG signal is local groups of pyramidal neurons of the cortex firing synchronously and is arranged perpendicular to its surface [60]. In our case, a subset of the lead field matrix was used for the 64 selected EEG electrodes.

Table 3: Summary presentation and description of the most important models and connectivity metrics or measures that were used in the methodological section of this study.

Metric or model name
Acronym
Description

Current cortical density
CCD
Forward problem solution
A model that aims to explain the correspondence of cortical source activity to scalp electrical potentials, taking account of skull and scalp conductivity.

Weighted minimum norm estimate
wMNE
Inverse problem solution
An estimation of how signals captured at the scalp correspond to source activations, with their power limited by the cortical physiology.

Directed transfer function
DTF
Granger causality measure
A metric of effective network connectivity (functional connectivity that incorporates causal relations instead of statistical inference alone) that produces directed networks with weighted edges.

Characteristic path length
CPL
Network integration
A representative measure of shortest distances between network nodes that are connected to each other.

Clustering coefficient
CC
Network segregation
A measure of the tendency of network nodes to become organized into functionally separated clusters.

Density
D
Network density
A measure of existing connections against the theoretical maximum number of possible connections if the network was fully connected.

Small-worldness
SW
Overall network effectiveness
A model of network behavior, where short paths and increased forming of functional clusters lead to optimization and resilience of information transfer.

Out-strength and in-strength
OS and IS

The total nodal sum of weights from outgoing and incoming connections, respectively.

Weighted minimum norm estimate was used to solve the ill-posed inverse problem (Table 3) by minimizing the source space energy based on the fact that the power of the source dipoles is limited by the cortex physiology [61]. Minimum norm estimate aims at minimizing the following equation:

where refers to the actual recordings from the scalp, is the simultaneous current dipoles to be calculated, is the lead field matrix, is the regularization parameter, and is the regularization term which in our case refers to the energy of the solution’s dipoles. The first term in the above equation represents the error between the actual and predicted electrode recordings. The second term is the penalization term, which aims at enforcing the abovementioned energy restriction. Lastly, , which balances the effect of the penalization term, was calculated using the L-curve method [62]. The solution of minimum norm estimate was derived using Tikhonov regularization in the regularization toolbox [63]. 24 custom-defined ROIs were created at the surface of the cortex model, in order to proceed to connectivity analysis, as illustrated in Figure 2. The ROI time series signal (Figure 3) was calculated by averaging the amplitude from all included cortical current dipoles.

2.2.3. Functional Connectivity

In total, 24 ROIs were defined on the cortical source model, consisting of the following areas bilaterally: presupplementary motor area (pSMA), supplementary motor area (SMA), dorsal premotor area (PMd), ventral premotor area (PMv), cingulate motor area (CMA), primary foot motor area (M1F), primary hand motor area (M1H), primary lip motor area (M1L), primary foot somatosensory area (S1F), primary hand somatosensory area (S1H), secondary somatosensory area (S2), and somatosensory association area (SAC) (Figure 2). Their average activation time series were computed for every time interval.

Directed transfer function (DTF) [64] was used (Table 3) in order to calculate functional cortical connectivity of the sensorimotor network consisting of the 24 ROIs as nodes [44], computing causal relations among the nodal activation time series. The produced connectivity matrices were thresholded, using the surrogate data method with testing of significance of connections, instead of using absolute or relative thresholding [65–67]. During computation of DTF, a number of 1000 surrogate permutations and a significance level of 0.05 were set, resulting in partially connected matrices with only the most significant causal connections. DTF is a measure based on Granger causality [68] that uses the multivariate autoregressive model described by the following function:

where is the model order, contains the ROIs values at time , is the residual noise vector, and is a coefficient k × k-sized matrix [44]. Using the above equation, the matrices are computed by means of the minimalization of the residual noise .

The order of the multivariate autoregressive model [69] was chosen to be 8 after considering the following criterions [70]: (a) tests that demonstrated an optimal order of 10 for a sampling rate of 128 Hz for modeling EEG spectra [71, 72]; (b) the model order should be smaller than  × Fs, where is the expected lag between two brain processes and Fs the sampling rate; (c) better to err on the side of selecting a larger model order, (d) using the same model order for all DTF computations.

Equation (3) is described in the frequency domain as

where is a transfer matrix of the system, and it contains information about the relationships between signals. It is nonsymmetric, so it allows for finding causal dependencies. DTF is then computed by the equation:

DTF describes casual influence of channel j on channel i at frequency f. For our analysis, DTF was computed for the networks formed at the frequency bands of alpha rhythm at 8–12 Hz (“alpha networks”) and beta rhythm at 13–30 Hz (“beta networks”), as those are considered the brainwaves most relevant to the sensorimotor processes [73].

2.2.4. Network Analysis

Network analysis in terms of descriptors of the weighted directed graphs (“network properties”) was performed with the brain connectivity toolbox for MATLAB [74] for the alpha and beta networks formed during early and late time intervals. Out-strength (OS), in-strength (IS), and clustering coefficient were computed for each of the 24 nodes of the network during each task of imagery. Characteristic path length (CPL), mean clustering coefficient (CC), and density (D) were also computed at the level of graphs of each task of imagery. Topology of small-worldness (SW) was then derived from these network properties for each task of imagery [44]. To facilitate further analysis, all properties were averaged for the 8 imagery categories mentioned in Experimental Procedure. A summary of the interpretation of these network properties is also presented in Table 3.

CPL, a measure of integration, calculates the sum of the shortest distances among connected graph nodes, divided by the number of nodes [75, 76]. It is described by (6), where is the average distance between node and the other node, and is the shortest path between nodes and . The distance matrix was computed using the logarithmic conversion of weights by the Floyd-Warshall algorithm [77], as implemented in the BCT [78].

Global CC, a measure of segregation, estimates the tendency of the graph nodes to organize into clusters [79, 80]. It is described by (7), where is the clustering value for a node , is the node’s degree, and is the number of neighboring nodes that connect to each other in triangles around node .

D is the ratio of a network’s actual connections to the maximum possible connections. In our example, the graphs were only partially connected, and the connection weights were ignored, since all connections were considered significant as computed by DTF with testing for statistical significance () by surrogate data. It is therefore described by (8), where is the ensemble of the network’s connections, and is the ensemble of the network’s nodes.

SW is defined as the combination of short paths and high clustering in a network, when compared to random networks with comparable paths constructed by the same number of nodes and connections. This property has commanded attention as an important brain network characteristic that models the brain’s effective communication patterns [81–83]. It is described by (9), and in our study, the comparison of CPL and CC was made against 10,000 random networks. These random networks were directed graphs, with the same number of nodes and edges as the original, and they were generated using the brain connectivity toolbox [84]. CPL and CC were computed for each random network and compared against the original, in a process that was iterated 10,000 times to produce a range of SW values. The range of values was then averaged to produce a single robust value of the SW property for each original network.

The strength of node strength is the sum of the weights of connections to or from that node. The IS, therefore, is the sum of incoming connection weights, and the OS is the sum of outgoing connection weights [74]. They are described by (10) and (11), respectively, where is the weighted directed connectivity matrix, with a direction of , is the number of columns, and is the number of rows.

2.3. Statistical Analysis

Adjacency matrices computed by directed transfer function were compared between healthy and patient groups using the network-based statistic toolbox for MATLAB [85]. We performed the statistical analysis using the false discovery rate on the general linear model with t-test [86], a significance level of 0.05 and 50,000 permutations. Using a between-group design, we compared alpha and beta networks of healthy subjects to alpha and beta networks of SCI subjects for each imagery category, elaborating the comparisons for the effect of early and late intervals. We also compared alpha to beta networks of each category and time interval using a within-group design. Differences in networks were visualized using the BrainNet Viewer for MATLAB [87].

Statistical analysis of age, imagery capacity, and computed network properties was performed in IBM SPSS Statistics (version 23), and we set a significance level of 0.05 for all statistical tests. All variables were explored for normality assumption (healthy and SCI groups as grouping factor) using visual inspection of histograms, normal Q-Q plots and boxplots, skewness and kurtosis [88–90], and normality tests (Shapiro-Wilk test and Kolmogorov-Smirnov test) [91, 92]. Depending on normality assumption, different analyses were performed (paired t-tests or Wilcoxon signed-rank tests). Normality assumption was met for the variable age and for the VVIQ score for both groups. Independent sample t-tests were performed to reveal significant age and VVIQ differences between the two groups. As the groups did not differ either for age distributions (healthy-skewness: 0.407 (SE = 0.687), kurtosis: −1.418 (SE = 1.334); SCI-skewness: 0.651 (SE = 0.687), kurtosis: −0.752 (SE = 1.334)) or for their reported imagery capacity (VVIQ: , , and ), the rest of the statistical analysis was planned accordingly.

We planned within-group comparisons of brain network properties using as grouping factor the rhythm (alpha, beta). Differences of variables between the two rhythms were calculated for the categories of visual motor imagery tasks, separately at early and late time intervals. Subsequently, we calculated within-group comparisons of brain network properties using as grouping factor the time interval (early, late). Between-group comparisons were performed using the calculated differences of variables at the two time intervals with either independent samples t-tests or Mann–Whitney U tests.

Nodal strengths, both incoming and outgoing, were averaged across different motor imagery tasks, rhythms, and time intervals, and total nodal strengths were calculated. They were tested for normality assumption for both groups and analyzed within groups using descriptive statistics and between groups (SCI, healthy) using Mann–Whitney U test. Targeted differences between nodes CMA_L and CMA_R were tested for statistical significance using either Pearson or Spearman correlation coefficient depending on normality assumption.

3. Results

3.1. Functional Connectivity

Visualizations of connectivity maps for the two groups (SCI, healthy) were made using the eConnectome toolbox for alpha and beta networks during both time intervals averaged across the motor imagery categories (Figure 4). The highest information transfer in all examined networks came from the bilateral cingulate motor areas including their reciprocal communication. In both groups, the maximum incoming nodal strength was observed in right CMA (Table 4), whereas the maximum outgoing nodal strength was found in left CMA (Table 5). Between-group comparisons revealed significant differences in total nodal strengths, both incoming and outgoing.

Figure 4: Average information transfer (calculated by directed transfer function) of healthy and SCI groups calculated for alpha (left) and beta (right) rhythm networks during the early imagery time interval, for the hands motor imagery category. Connections between the bilateral cingulate motor areas (CMA_R ← → CMA_L) presented the highest information transfer. Only connections with at least 25% of max information transfer among all statistically significant connections are displayed.

Table 4: Descriptive statistics of nodal incoming strengths with statistically significant differences between healthy and SCI participants and between-group comparison results.

Nodes
Median
Mean ranks
IQR: [Q₁, Q₂]
Healthy versus SCI

Healthy
SCI
Healthy
SCI
Healthy
SCI

SIF L
0.235
0.230
37.72
27.28
[0.231, 0.253]
[0.226, 0.238]
,

SIH L
0.263
0.246
44.13
20.88
[0.258, 0.270]
[0.239, 0.256]
,

SIH R
0.256
0.233
40.94
24.06
[0.243, 0.265]
[0.229, 0.248]
,

CMA L
0.410
0.438
18.84
46.16
[0.402, 0.419]
[0.431, 0.445]
,

CMA R
0.486
0.473
38.38
26.63
[0.475, 0.509]
[0.458, 0.489]
,

SMA L
0.236
0.217
41.34
23.66
[0.230, 0.243]
[0.213, 0.229]
,

SMA R
0.305
0.289
39.50
25.50
[0.296, 0.314]
[0.281, 0.297]
,

PMd L
0.309
0.272
46.84
18.16
[0.296, 0.320]
[0.260, 0.278]
,

PMd R
0.299
0.284
37.78
27.22
[0.290, 0.304]
[0.274, 0.299]
,

PMv L
0.252
0.232
40.53
24.47
[0.236, 0.259]
[0.220, 0.243]
,

PMv R
0.262
0.247
41.81
23.19
[0.254, 0.267]
[0.240, 0.255]
,

Table 5: Descriptive statistics of nodal outgoing strengths with statistically significant differences between healthy and SCI participants and between-group comparison results.

Nodes
Median
Mean ranks
IQR: [Q₁, Q₂]
Healthy versus SCI

Healthy
SCI
Healthy
SCI
Healthy
SCI

S1F L
0.073
0.084
21.34
43.66
[0.080, 0.095]
[0.080, 0.095]
,

S1F R
0.131
0.133
31.81
33.19
[0.123, 0.146]
[0.123, 0.146]
,

S1H L
0.090
0.084
38.38
26.63
[0.084, 0.100]
[0.077, 0.091]
,

S1H R
0.091
0.166
16.50
48.50
[0.086, 0.096]
[0.153, 0.181]
,

SAC L
0.071
0.087
18.28
46.72
[0.069, 0.074]
[0.084, 0.092]
,

SAC R
0.025
0.015
47.41
17.59
[0.023, 0.027]
[0.014, 0.016]
,

S2 L
0.010
0.022
17.41
47.59
[0.009, 0.012]
[0.020, 0.024]
,

S2 R
0.128
0.148
20.91
44.09
[0.124, 0.138]
[0.139, 0.159]
,

M1F L
0.078
0.073
38.13
26.88
[0.072, 0.088]
[0.069, 0.076]
,

M1F R
0.294
0.236
47.69
17.31
[0.286, 0.324]
[0.230, 0.263]
,

M1H L
0.311
0.170
47.59
17.41
[0.283, 0.339]
[0.163, 0.181]
,

M1H R
0.061
0.059
38.34
26.66
[0.059, 0.066]
[0.058, 0.061]
,

M1L L
0.003
0.003
26.25
38.75
[0.003, 0.003]
[0.003, 0.004]
,

M1L R
0.011
0.010
42.78
22.22
[0.010, 0.013]
[0.009, 0.010]
,

CMA L
2.927
2.519
47.28
17.72
[2.805, 3.008]
[2.404, 2.581]
,

CMA R
0.879
1.156
18.41
46.59
[0.784, 0.973]
[1.082, 1.261]
,

SMA L
0.350
0.276
46.25
18.75
[0.326, 0.375]
[0.267, 0.295]
,

SMA R
0.206
0.187
41.72
23.28
[0.191, 0.226]
[0.176, 0.194]
,

pSMA L
0.044
0.031
46.50
18.50
[0.040, 0.048]
[0.030, 0.035]
,

pSMA R
0.026
0.031
19.28
45.72
[0.024, 0.027]
[0.030, 0.033]
,

PMd L
0.149
0.187
21.13
43.88
[0.141, 0.167]
[0.179, 0.201]
,

PMd R
0.484
0.577
18.25
46.75
[0.448, 0.521]
[0.559, 0.606]
,

PMv L
0.014
0.024
16.50
48.50
[0.013, 0.014]
[0.022, 0.025]
,

PMv R
0.016
0.032
17.91
47.09
[0.014, 0.018]
[0.028, 0.035]
,

More precisely, significant differences in incoming strengths were found bilaterally in CMA, SMA, S1H, PMd, and PMv as well as in the left S1F. Healthy participants showed higher incoming strengths in all aforementioned nodes apart from the left CMA (Figure 5) compared to the SCI participants (Table 4).

(a)

(b)

Figure 5: Nodal strengths (a: out-strength, b: in-strength) of bilateral cingulate motor areas for both subject groups. Left cingulate motor area (CMA_L) showed the highest out-strength and right cingulate motor area (CMA_R) showed the highest in-strength in the network. SCI subjects presented significantly higher CMA_R out-strength and CMA_L in-strength than healthy subjects. “” represent extreme values and “o” represent outliers.

Outgoing strengths were found to be significantly different between groups in all nodes (Table 5) apart from the right S1F. In more detail, SCI group showed considerably higher outgoing strengths in the S1F and SAC in the left hemisphere, in S1H and CMA in the right hemisphere as well as in S2, PMd and PMv bilaterally. In the remaining nodes, healthy participants were found to have increased outgoing strengths compared to SCI participants.

When calculating differences of group averages (healthy group-SCI group) of in-strength (IS) and out-strength (OS) of cingulate motor areas (CMAs) during all motor imagery categories, a possible trend was revealed. OS of CMA_R was consistently higher in the healthy group, while OS of CMA_L was consistently higher in the SCI group. The opposite held true for IS of those nodes (Figure 6). Between-group differences in CMA_R were negatively correlated to those in CMA_L in targeted imagery categories (early alpha walking (, ), late alpha walking (), early beta walking (, ), and late beta walking (, )).

Figure 6: Differences of group averages (healthy-SCI) of in-strength (IS) and out-strength (OS) of cingulate motor areas (CMAs) during all motor imagery categories. A trend was revealed, in which OS of CMA_R was consistently higher in the healthy than the SCI group. OS of CMA_L was consistently lower in the healthy than the SCI group. The opposite held true for IS of those nodes. This trend reached statistical significance only for early alpha walking () and late beta walking () tasks.

3.2. Network-Based Statistics between Groups and within Groups

Important differences of connectivity were found only in between groups (SCI against healthy), where a subset of connections had significantly higher FC in the healthy group than in the SCI group in the hands motor imagery category (Figure 7). This subset included connections with lower FC in the SCI group of M1H_R to bilateral primary foot motor areas (M1F), primary foot and hand sensory areas (S1H, S1F), the somatosensory association areas (SAC), and the secondary sensory areas (S2). This finding persisted in both alpha and beta networks when testing with t-test.

Figure 7: From the comparison of the networks of healthy and SCI subjects, a subset of network connections emerged as significantly stronger in the healthy group than in the SCI group for both the alpha and beta networks of “hands” motor imagery category, as calculated by network-based statistics—false discovery rate methodology.

When we further tested the networks of our participants by grouping the SCI subjects by outcome (positive and negative), no differences were found between the networks of the two groups. During all permutations, the value of the false discovery rate did not approach statistical significance (). Also, testing for other imagery categories did not reveal significant differences, with the exception of the walking category. Comparing the networks of healthy and patients during the walking imagery category, significantly greater S2_L-PMv_R connectivity was found in the SCI group.

Furthermore, when testing for main effect of within time interval and brainwave rhythm within the healthy and patient groups, no further statistical significant differences of the connectivity weights of the network were observed.

3.3. Analysis of Network Properties

3.3.1. Within-Group Comparisons of Graph Properties between Alpha and Beta Showed Less Segregation, Less Integration, Greater Density, and Less Effectiveness of Beta Networks

When exploring within-group differences of graph properties using as grouping factor the rhythm (alpha, beta), beta networks showed less segregation, less integration, and less overall effectiveness compared to alpha networks. CPL, CC, and SW showed significantly lower values in beta compared to alpha networks, in both early and late time intervals. These findings were observed during nearly all imagery categories in both SCI and healthy group.

On the opposite, beta networks showed greater density compared to alpha networks. D was significantly greater in beta networks in both early and late time intervals of all imagery categories in both SCI and healthy group (Figure 8). Aggregated statistical test results and values for abovementioned findings can be found in supplementary material (available here). Specific exemptions are detailed below, as differences of graph properties in beta network compared to alpha did not reach statistical significance only in walking category, in the following cases: (a) in SCI group, CPL during the late interval, (b) in SCI group, CC and SW during the early interval, and (c) in healthy group, CC during the late interval.

Figure 8: Network density was significantly higher in beta than alpha networks (all values < 0.001) of the “hands” motor imagery category during both early and late time intervals for all subjects of both the SCI and the healthy group (all values < 0.001). This finding was also consistent across every other motor imagery categories in both groups (all values < 0.001).

3.3.2. Within-Group Comparisons of Graph Properties between Early and Late Time Intervals Showed (1) Less SCI Network Integration during Late Walking Imagery, (2) Greater SCI Network Segregation and Stable Effectiveness for Distal Tasks during Late Imagery, and (3) Less Healthy Network Segregation and Effectiveness for Distal Tasks during Late Imagery

When exploring within-group differences of graph properties using as grouping factor the time interval (early, late), few significant differences were observed.

Regarding network integration, significant difference of CPL values was shown only SCI group’s alpha networks during the walking task (, , and ). More precisely, SCI subjects were characterized by lower path lengths at the second stage of the task (late) (early alpha walking CPL: 7.809; late alpha walking CPL: 7.032). Changes in the CPL in the beta rhythm comparing the two time intervals were not observed. Also, no significant difference was found for healthy subjects.

Regarding network segregation, significantly higher CC values in the SCI group were observed in the alpha band of distal imagery category (, , and ; early alpha distal CC: 0.0076; late alpha distal CC: 0.0082). Considerable differences in mean CC at the beta band were not found. In healthy participants, considerably lower CC value was found only in alpha band of the left category (, , and ; early alpha left CC: 0.0094; late alpha left CC: 0.0086).

Regarding network density, healthy group showed significantly higher D values at the late stage of linear imagery tasks in alpha rhythm (, , and ; early alpha linear D: 0.3595; late alpha linear D: 0.3713), whereas density was considerably less at the late stage of proximal tasks in beta rhythm (, , and ; early beta proximal D: 0.5904; late beta proximal D: 0.5784). SCI group showed greater density when comparing the two time intervals of right imagery category in alpha band (, , and ; early alpha right D: 0.3663; late alpha right D: 0.3801), but no alterations were found in beta networks.

Regarding overall network effectiveness, significant results of SW were found for healthy subjects at distal imagery tasks in both alpha and beta rhythms (alpha: , , and ; beta: , , and ). In more detail, significantly lower values of SW were observed in distal imagery tasks between the two time intervals (early alpha distal SW: 1.553; late alpha distal SW: 1.406; early beta distal SW: 1.159; late beta distal SW: 1.137). For the SCI group, considerable differences were not observed.

3.3.3. Between-Group Comparisons of Graph Properties Showed Not Only Similar Network Integration and Density But Also Greater Segregation and Effectiveness of Alpha Band Networks in Some Imagery Categories for the Patients

When exploring between-group differences of graph properties, few significant differences were observed. Comparisons of CPL and D did not reveal any considerable difference across any imagery category (supplementary material), showing similar network integration and network density of the networks of healthy and patient subjects.

Regarding network segregation, significant changes in CC were observed at the left imagery tasks (, , and ), at the rotational imagery tasks (, , and ), and at distal imagery tasks (, ), all appearing in alpha band. In more detail, SCI subjects seem to show greater CC of alpha networks in the late part of the aforementioned imagery tasks (alpha left dif (mean)—SCI: 0.00019, healthy: −0.00084; alpha distal dif (median)—SCI: 0.00049, healthy: −0.00075; alpha rotational dif (mean)—SCI: 0.00061, healthy: −0.0010). Significant differences were not found between groups in the beta networks of all tasks.

Regarding overall network effectiveness, the SCI group seems to have only a significant change in mean SW of alpha network while performing VMI on distal imagery tasks compared to healthy (, , and ). More precisely, SCI group seems to show greater SW of alpha networks in the late part of the distal tasks (alpha distal dif (mean)—SCI: 0.0253, healthy: −0.1467). Other between-group differences were not observed (supplementary material).

4. Discussion

4.1. Functional Connectivity

A subset of the sensorimotor network during hands motor imagery was shown to have significantly lower functional connectivity power in the SCI group compared to the healthy group, a finding from analyses of the general linear model. This subset included connections of the M1H_R cortical area (theoretically the nondominant hand primary motor area for right-handed subjects) with other motor and sensory cortical areas. This subset also excluded the “assistive” motor nodes (CMA/SMA/pSMA/PMv/PMd). Interestingly, among these excluded nodes were also the ones that were shown to have consistently higher OS and IS, for all imagery categories, as we will discuss later on. The subset of connected nodes to M1H_R included bilateral primary foot motor areas (M1F), primary foot and hand sensory areas (S1H, S1F), the somatosensory association areas (SAC) located in superior parietal cortex (SPC), and the secondary sensory areas (S2). Small differences between alpha and beta rhythm networks can be observed, whereas most of those connections’ lower FC reached statistical significance for either rhythm. These cortical areas can be identified as the point of origin of the pyramidal tract and the point of conclusion of the major somatosensory tracts. This finding could suggest that chronic disruption of reciprocal communication between the brain and spinal cord, even in noncomplete injuries, could result in permanent significant decrease of connectivity between a subset of the functional sensorimotor network at the cortical level. This effect was observed regardless of positive or negative neurological outcome since grouping SCI subjects by outcome did not reveal any differences regarding this finding. While the lack of difference between those two clinically and functionally different subgroups of patients could be affected by a lack of statistical power when comparing small samples, a possible explanation could also be explored along the lines of mental imagery capacity. Motor execution and motor imagery do present differences in the level of neural networks that are affected by the subjects’ quality and intensity of imagination [93]. In our study, the patients did not report differences in vividness of visual imagery to that of the healthy participants, as measured by the VVIQ questionnaire at the time of the experimental procedure. A degree of “though extinction process” has been reported in chronic paralysis [94], but since this is not the case in our investigation, it is possible that the lack of difference in the networks of patients with positive and negative outcomes could be due to the unaffected mental capacity of the participants.

Further significant differences between the healthy and patient groups were found with regard to nodal strengths. Almost all nodes had significantly different out-strength between the two groups. An interesting pattern can be identified. Primary motor areas, supplementary motor areas, left presupplementary motor area, and left cingulate motor area show significantly higher out-strength in the healthy group. On the contrary, premotor areas, right presupplementary motor area, and right cingulate motor area show significantly higher out-strength in the SCI group. Incoming connection strengths to primary somatosensory areas, premotor areas, supplementary motor areas, and right cingulate motor area are higher for healthy participants, while in-strength of left cingulate motor area is higher for SCI participants. Despite the differences in strengths of those nodes, similar patterns of connectivity were found for both groups [12, 13]. The significantly reduced nodal strengths could reflect the disconnection itself and the reduced input and output of the sensorimotor pathways in spinal cord injury. Nonetheless, the higher in-strengths of left cingulate motor area and out-strengths of premotor areas and right pSMA and CMA in SCI group could indicate an attempt of the sensorimotor network to compensate for the impaired function [10, 11, 95].

CMA areas have been previously identified as important information hubs for sensorimotor networks, especially those of beta rhythm [9, 11]. In our study, this attribute is confirmed, since bilateral CMA areas consistently received the greater inflow and contributed the greater outflow in terms of connection strengths for all categories of motor imagery. Moreover, their reciprocal communication constituted the most powerful connections of every examined network. On the other hand, an important role has been identified for the SMAs [11–13] that have been shown to present notable outflow during motor imagery tasks and form clusters with the CMAs. Their role was asserted in our work previously too [44, 96], but it is not so apparent in our current study, where the SMAs were not among the top contributors in either outflow or inflow. Although not easily explained, the meaning of this finding can be explored along possible factors: (a) the random-oddball (unexpected imagery task) paradigm of presentation, (b) the MVAR model order set, and (c) the definitions of the midline network nodes themselves. The degree that each factor possibly contributed to this finding is an issue for further investigation. An example of SMAs and primary motor areas not presenting the greatest strength during hands motor imagery has also been recently reported in a study [97] where the authors used transcranial direct current stimulation to affect the connectivity of a broader definition of sensorimotor ROIs.

With regard to differentiating different upper limb motor imagery tasks, our results did not produce significant differences in terms of spatial patterns specific to certain tasks. Moreover, network-based differences between healthy and patients, although significant for the all-inclusive upper limb imagery category, did not reach statistical significance for specific categories, suggesting possibly a lack of statistical power for these categorical differences. This is not unexpected, since connectivity features, in general, have so far shown only moderate success in classification of motor imagery tasks [98, 99]. It should be also noted that some effort has been made in analyzing effective networks of compound motor imagery tasks [100]. Differentiating anatomical levels and consecutive classification should perhaps be better explored along the lines of time-varying connectivity [95, 101–103] instead of spatial pattern analysis.

Walking motor imagery, while it also did not reveal specific connectivity patterns, produced the most promising results in terms of classification, in accordance to previous studies suggesting that maximally different conditions should be explored [98]. Walking motor imagery category was the only one where the negative correlation between-group differences of the two cingulate motor area strengths reach statistical significance in half of the studied cases, those of early alpha rhythm walking networks () and of late beta walking networks (). Moreover, the comparison of networks of healthy and patient subjects produced at least one significantly stronger connection, between the right ventral premotor area and the left secondary somatosensory area, although it is unclear whether this can be attributed to plasticity or merely to SCI-induced disconnection sequelae. To the best of the authors’ knowledge, this is the first electroencephalographic study of functional cortical connectivity after incomplete spinal cord injury, and it is also the first functional cortical connectivity study examining multiple motor imagery tasks in those patients regardless of recording modality.

4.2. Analysis of Network Properties

Analysis of within-group effect of rhythm produced the most consistent results. According to the revealed pattern, alpha networks present lower integration (as measured by CPL), higher segregation (as measured by CC), while being less dense and more “effective” (as measured by SW) than beta networks. These findings are present across all motor imagery categories and they closely match findings from our previous study on the role of alpha and beta rhythms in sensorimotor networks [44, 104]. Our previous work suggested a pattern where alpha rhythm engaged local information processing using greater wiring costs [105], and that beta rhythm assumes a coordinative role during the sensorimotor process [106]. These findings were then observed on different ROI models and during simpler but far more repetitive motor imagery and motor execution tasks. They are also now replicated on a wider definition of the model of ROIs and during multiple, more complex, random motor imagery tasks. More importantly, these findings have now also been confirmed on networks of SCI patients with incomplete injury, allowing us to attempt to model the behavior of other between-group and within-group findings based on this pattern of alpha and beta organization.

Between the two groups, the fact that CPL and D were not significantly different, neither in alpha nor in beta networks, allows us to make direct comparisons of their sensorimotor network organization since they reach the same level of wiring costs and node integration. Moreover, the few between-group differences were observed mostly in alpha rhythm, which could be interpreted as differences only in local processing in the sensorimotor network of SCI patients. Increased functional segregation (CC in left, rotational and distal categories) and increased “effectiveness” of the network (SW in distal) were found for certain categories of motor imagery. More importantly, they were observed for distal arm imagery tasks, those that correspond to spinal cord levels below the level of injury, as the majority of the SCI subjects included in our study suffered from mid to low cervical SCI (C4–C8). As those differences were also observed during the late time interval, they can possibly be interpreted as an effect for delayed adaptation (compensation) of the sensorimotor network at the cortical level. This could possibly fall in line with reported increased network fault tolerance [10] and an increase of local efficiency and communication between closest cortical areas [15] during paralyzed foot motor imagery that has been reported in chronic complete SCI subjects.

Regarding the walking imagery category in our study, walking networks were the only where alpha and beta rhythm differences did not reach statistical significance in certain cases. Indeed, this previously reported increase of local efficiency and close communication in complete SCI appears to also be possible in incomplete injuries as well. The walking beta networks did not show less segregation and effectiveness than the walking alpha networks of incomplete SCI subjects during the early imagery part and also did not show less integration during the late imagery part, as was the case with the upper limb imagery categories. These findings suggest the presence of a phenomenon that has been attributed to adaptive plasticity and compensation when it regards patients with complete injuries [10, 15]. Within-group effect of time interval for upper limb motor imagery categories was far more sporadic, showing greater network segregation and less network integration of the alpha rhythm network in the patient group during the late imagery in some categories. What could interestingly fall into place with the rest of the interpretation is that healthy subjects display a drop of overall network effectiveness (as depicted by SW) in both alpha and beta networks during the late imagery part. In the SCI group, such a difference was not observed, an observation that could possibly be attributed to the same delayed compensatory effect induced by the injury. Since indeed the rest of the observed effects are not consistent, our reported findings cannot be obligatory attributed to neuroplasticity effects. Therefore, it is evident that more investigation in the direction of modeling the effect of spinal cord injury on the effective connectivity of the brain during different time points and injury severity conditions is needed before drawing accurate conclusions.

4.3. Limitations and Future Work

Among limitations of our study, one investigating functional connectivity should remain wary when the analysis reveals significant differences between the groups, as those differences are usually very small and not always clear if functionally or clinically meaningful. As such, it remains difficult to identify compelling advantages of graph theory-based analysis of brain activity over other approaches to provide additional important insight into the effects of SCI on brain activity. Functional connectivity at the source level also suffers from certain disadvantages including localization error, smoothing effect, and a degree of uncertainty of the connectivity between spatially close nodes [107]. There are several factors contributing to source localization errors, induced most importantly by the forward model but also by the inverse. The resolution of the source reconstruction is determined by the source space, with 4.69 mm of average source spacing and 22.04 mm² average surface area per source. Nonetheless, resolution of the source space is considered sufficient for the purpose of the study, and the same model has been previously used by similar studies [9–15]. The 3-layer block element modifier forward model introduces errors through the use of a general template anatomy, modeling of skull conductivity as isotropic, not modeling cerebrospinal fluid and conductivity ratios. Moreover, it is known that EEG boasts great temporal but suffers from low spatial resolution and has been traditionally considered able to detect rapid brain dynamics in a trade-off with source estimation and low signal to noise ratio due to volume conduction effect [108, 109]. EEG in general greatly suffers from anisotropic conductivity of skull leading to signals blurring and to low EEG spatial resolution. In accordance, the localization error of deeper sources is considered greater than swallower ones. An interesting approach to address these problems would be the investigation of time-adaptive connectivity with a focus on temporal alterations of important connections rather than spatial [64], while using individual subject anatomy, an approach that we aim to explore in our future work.

5. Conclusions

We observed that chronic disruption of reciprocal communication between the brain and spinal cord, even in the context of incomplete injuries, could result in permanent significant decrease of connectivity between a subset of the functional sensorimotor network at the cortical level. This effect was observed regardless of positive or negative neurological outcome since grouping SCI subjects by outcome did not reveal any further difference. Cingulate motor areas were identified as important information hubs in different categories of motor imagery as they consistently showed the highest in-strengths (CMA_L) and out-strengths (CMA_R) in both groups of participants. While SCI subjects also followed the same pattern, they had higher outflow from left CMA and higher inflow to right CMA than healthy subjects. For both groups, alpha networks were less dense while having both longer average paths and more clustering than beta networks in almost all imagery categories. SCI patients showed signs of increased local processing in the late part of imagery, possibly an adaptive compensatory mechanism of injury-induced neuroplasticity.

Ethical Approval

This study was conducted in accordance with the Declaration of Helsinki (1964) and its following amendments. The Bioethics & Ethics Committee of Faculty of Medicine, Aristotle University of Thessaloniki, approved the study.

Consent

All experiments were conducted with the subjects’ understanding and written informed consent.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this article.

Acknowledgments

Early phases of the study were supported in part by a Mario Boni grant from the Cervical Spine Research Society-European Section, awarded to CSI: Brainwave Project. This study was supported by 2018 “Brihaye” EANS Research Fund Grant from the European Association of Neurosurgical Societies. This work was supported by European Union’s Horizon 2020 research and innovation programme under grant agreement No 681120 of the SmokeFreeBrain project (http://smokefreebrain.eu).

Supplementary Materials

Statistical analysis: methodology and results for CPL, CC. D, SW, and correlations. Figure 1: SCI group shows increase in mean global clustering coefficient compared to healthy group in late time interval of the alpha networks, during “left,” “distal,” and “rotational” motor imagery categories. Figure 2: SCI group shows increased small-worldness (SW) compared to healthy group in the late time interval of alpha networks for “distal” imagery category. Figure 3: between-group differences in CMA_L were negatively correlated to those in CMA_R in early alpha walking and late beta walking. Table 1: descriptive statistics of total nodal out-strengths and in-strengths of the networks of the SCI subjects. Table 2: descriptive statistics of total nodal out-strengths and in-strengths of the networks of the healthy subjects. (Supplementary Materials)

References

K. D. Anderson, M. E. Acuff, B. G. Arp et al., “United States (US) multi-center study to assess the validity and reliability of the spinal cord independence measure (SCIM III),” Spinal Cord, vol. 49, no. 8, pp. 880–885, 2011.

R. Nardone, Y. Höller, F. Brigo et al., “Descending motor pathways and cortical physiology after spinal cord injury assessed by transcranial magnetic stimulation: a systematic review,” Brain Research, vol. 1619, pp. 139–154, 2015.

A. Athanasiou, M. A. Klados, N. Pandria et al., “A systematic review of investigations into functional brain connectivity following spinal cord injury,” Frontiers in Human Neuroscience, vol. 11, p. 517, 2017.

R. Nardone, Y. Höller, F. Brigo et al., “Functional brain reorganization after spinal cord injury: systematic review of animal and human studies,” Brain Research, vol. 1504, pp. 58–73, 2013.

P. Freund, N. Weiskopf, J. Ashburner et al., “MRI investigation of the sensorimotor cortex and the corticospinal tract after acute spinal cord injury: a prospective longitudinal study,” The Lancet Neurology, vol. 12, no. 9, pp. 873–881, 2013.

J. Hou, Z. Xiang, R. Yan et al., “Motor recovery at 6 months after admission is related to structural and functional reorganization of the spine and brain in patients with spinal cord injury,” Human Brain Mapping, vol. 37, no. 6, pp. 2195–2209, 2016.

W. Zheng, Q. Chen, X. Chen et al., “Brain white matter impairment in patients with spinal cord injury,” Neural Plasticity, vol. 2017, Article ID 4671607, 8 pages, 2017.

T. Ilvesmäki, E. Koskinen, A. Brander, T. Luoto, J. Öhman, and H. Eskola, “Spinal cord injury induces widespread chronic changes in cerebral white matter,” Human Brain Mapping, vol. 38, no. 7, pp. 3637–3647, 2017.

L. Astolfi, H. Bakardjian, F. Cincotti et al., “Estimate of causality between independent cortical spatial patterns during movement volition in spinal cord injured patients,” Brain Topography, vol. 19, no. 3, pp. 107–123, 2007.

F. D. V. Fallani, L. Astolfi, F. Cincotti et al., “Cortical functional connectivity networks in normal and spinal cord injured patients: evaluation by graph analysis,” Human Brain Mapping, vol. 28, no. 12, pp. 1334–1346, 2007.

F. De Vico Fallani, L. Astolfi, F. Cincotti et al., “Extracting information from cortical connectivity patterns estimated from high resolution EEG recordings: a theoretical graph approach,” Brain Topography, vol. 19, no. 3, pp. 125–136, 2007.

D. Mattia, F. Cincotti, L. Astolfi et al., “Motor cortical responsiveness to attempted movements in tetraplegia: evidence from neuroelectrical imaging,” Clinical Neurophysiology, vol. 120, no. 1, pp. 181–189, 2009.

R. Sinatra, F. de Vico Fallani, L. Astolfi et al., “Cluster structure of functional networks estimated from high-resolution EEG data,” International Journal of Bifurcation and Chaos, vol. 19, no. 2, pp. 665–676, 2009.

L. Astolfi, F. Cincotti, D. Mattia et al., “Time-varying cortical connectivity estimation from noninvasive, high-resolution EEG recordings,” Journal of Psychophysiology, vol. 24, no. 2, pp. 83–90, 2010.

F. De Vico Fallani, F. A. Rodrigues, L. da Fontoura Costa et al., “Multiple pathways analysis of brain functional networks from EEG signals: an application to real data,” Brain Topography, vol. 23, no. 4, pp. 344–354, 2011.

J.-M. Hou, T.-S. Sun, Z.-M. Xiang et al., “Alterations of resting-state regional and network-level neural function after acute spinal cord injury,” Neuroscience, vol. 277, pp. 446–454, 2014.

Y.-S. Min, J. W. Park, S. U. Jin et al., “Alteration of resting-state brain sensorimotor connectivity following spinal cord injury: a resting-state functional magnetic resonance imaging study,” Journal of Neurotrauma, vol. 32, no. 18, pp. 1422–1427, 2015.

Y.-S. Min, Y. Chang, J. W. Park et al., “Change of brain functional connectivity in patients with spinal cord injury: graph theory based approach,” Annals of Rehabilitation Medicine, vol. 39, no. 3, pp. 374–383, 2015.

A. Oni-Orisan, M. Kaushal, W. Li et al., “Alterations in cortical sensorimotor connectivity following complete cervical spinal cord injury: a prospective resting-state fMRI study,” PLoS One, vol. 11, no. 3, article e0150351, 2016.

M. Kaushal, A. Oni-Orisan, G. Chen et al., “Evaluation of whole-brain resting-state functional connectivity in spinal cord injury: a large-scale network analysis using network-based statistic,” Journal of Neurotrauma, vol. 34, no. 6, pp. 1278–1282, 2017.

M. Kaushal, A. Oni-Orisan, G. Chen et al., “Large-scale network analysis of whole-brain resting-state functional connectivity in spinal cord injury: a comparative study,” Brain Connectivity, vol. 7, no. 7, pp. 413–423, 2017.

A. H. Hawasli, J. Rutlin, J. L. Roland et al., “Spinal cord injury disrupts resting-state networks in the human brain,” Journal of Neurotrauma, vol. 35, no. 6, pp. 864–873, 2018.

H. Luo, W. Dou, Y. Pan et al., “Joint analysis of multi-level functional brain networks,” in 2016 9th International Congress on Image and Signal Processing, BioMedical Engineering and Informatics (CISP-BMEI), pp. 1521–1526, Datong, China, 2016.

Y. Pan, W. Dou, Y. Wang et al., “Non-concomitant cortical structural and functional alterations in sensorimotor areas following incomplete spinal cord injury,” Neural Regeneration Research, vol. 12, no. 12, pp. 2059–2066, 2017.

S. Hétu, M. Grégoire, A. Saimpont et al., “The neural network of motor imagery: an ALE meta-analysis,” Neuroscience & Biobehavioral Reviews, vol. 37, no. 5, pp. 930–949, 2013.

M. Lotze and U. Halsband, “Motor imagery,” Journal of Physiology-Paris, vol. 99, no. 4–6, pp. 386–395, 2006.

S. Kraeutner, A. Gionfriddo, T. Bardouille, and S. Boe, “Motor imagery-based brain activity parallels that of motor execution: evidence from magnetic source imaging of cortical oscillations,” Brain Research, vol. 1588, pp. 81–91, 2014.

L. Xu, H. Zhang, M. Hui et al., “Motor execution and motor imagery: a comparison of functional connectivity patterns based on graph theory,” Neuroscience, vol. 261, pp. 184–194, 2014.

A. Guillot, C. Collet, V. A. Nguyen, F. Malouin, C. Richards, and J. Doyon, “Brain activity during visual versus kinesthetic imagery: an fMRI study,” Human Brain Mapping, vol. 30, no. 7, pp. 2157–2172, 2009.

N. Mizuguchi, M. Nakamura, and K. Kanosue, “Task-dependent engagements of the primary visual cortex during kinesthetic and visual motor imagery,” Neuroscience Letters, vol. 636, pp. 108–112, 2017.

J. Harris and A. Hebert, “Utilization of motor imagery in upper limb rehabilitation: a systematic scoping review,” Clinical Rehabilitation, vol. 29, no. 11, pp. 1092–1107, 2015.

L. Bunketorp Käll, R. J. Cooper, J. Wangdell, J. Fridén, and M. Björnsdotter, “Adaptive motor cortex plasticity following grip reconstruction in individuals with tetraplegia,” Restorative Neurology and Neuroscience, vol. 36, no. 1, pp. 73–82, 2018.

S. W. Tung, C. Guan, K. K. Ang et al., “Motor imagery BCI for upper limb stroke rehabilitation: an evaluation of the EEG recordings using coherence analysis,” in 2013 35th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), pp. 261–264, Osaka, Japan, 2013.

A. R. C. Donati, S. Shokur, E. Morya et al., “Long-term training with a brain-machine interface-based gait protocol induces partial neurological recovery in paraplegic patients,” Scientific Reports, vol. 6, no. 1, article 30383, 2016.

V. Rajasekaran, E. López-Larraz, F. Trincado-Alonso et al., “Volition-adaptive control for gait training using wearable exoskeleton: preliminary tests with incomplete spinal cord injury individuals,” Journal of NeuroEngineering and Rehabilitation, vol. 15, no. 1, p. 4, 2018.

A. Athanasiou, I. Xygonakis, N. Pandria et al., “Towards rehabilitation robotics: off-the-shelf BCI control of anthropomorphic robotic arms,” BioMed Research International, vol. 2017, Article ID 5708937, 17 pages, 2017.

G. Arfaras, A. Athanasiou, P. Niki et al., “Visual versus kinesthetic motor imagery for BCI control of robotic arms (Mercury 2.0),” in 2017 IEEE 30th International Symposium on Computer-Based Medical Systems (CBMS), pp. 440–445, Thessaloniki, Greece, 2017.

A. Athanasiou, G. Arfaras, I. Xygonakis et al., “Commercial BCI control and functional brain networks in spinal cord injury: a proof-of-concept,” in 2017 IEEE 30th International Symposium on Computer-Based Medical Systems (CBMS), pp. 262–267, Thessaloniki, Greece, 2017.

A. Athanasiou, G. Arfaras, N. Pandria et al., “Wireless brain-robot interface: user perception and performance assessment of spinal cord injury patients,” Wireless Communications and Mobile Computing, vol. 2017, Article ID 2986423, 16 pages, 2017.

“Brainwave control of a wearable robotic arm for rehabilitation and neurophysiological study in cervical spine injury (CSI:Brainwave),” https://clinicaltrials.gov/ct2/show/NCT02443558.

D. F. Marks, “New directions in mental imagery research,” Journal of Mental Imagery, vol. 19, no. 3-4, pp. 153–167, 1995.

R. Oostenveld and P. Praamstra, “The five percent electrode system for high-resolution EEG and ERP measurements,” Clinical Neurophysiology, vol. 112, no. 4, pp. 713–719, 2001.

O. Christ and M. Reiner, “Perspectives and possible applications of the rubber hand and virtual hand illusion in non-invasive rehabilitation: technological improvements and their consequences,” Neuroscience & Biobehavioral Reviews, vol. 44, pp. 33–44, 2014.

A. Athanasiou, M. A. Klados, C. Styliadis, N. Foroglou, K. Polyzoidis, and P. D. Bamidis, “Investigating the role of alpha and beta rhythms in functional motor networks,” Neuroscience, 2016, In press.

L. Qin, L. Ding, and B. He, “Motor imagery classification by means of source analysis for brain–computer interface applications,” Journal of Neural Engineering, vol. 1, no. 3, pp. 135–141, 2004.

B. Kamousi, A. N. Amini, and B. He, “Classification of motor imagery by means of cortical current density estimation and Von Neumann entropy,” Journal of Neural Engineering, vol. 4, no. 2, pp. 17–25, 2007.

B. Kamousi, Z. Liu, and B. He, “Classification of motor imagery tasks for brain-computer interface applications by means of two equivalent dipoles analysis,” IEEE Transactions on Neural Systems and Rehabilitation Engineering, vol. 13, no. 2, pp. 166–171, 2005.

B. J. Edelman, B. Baxter, and B. He, “EEG source imaging enhances the decoding of complex right-hand motor imagery tasks,” IEEE Transactions on Biomedical Engineering, vol. 63, no. 1, pp. 4–14, 2016.

R. Oostenveld, P. Fries, E. Maris, and J.-M. Schoffelen, “FieldTrip: open source software for advanced analysis of MEG, EEG, and invasive electrophysiological data,” Computational Intelligence and Neuroscience, vol. 2011, Article ID 156869, 9 pages, 2011.

A. Belouchrani and A. Cichocki, “Robust whitening procedure in blind source separation context,” Electronics Letters, vol. 36, no. 24, p. 2050, 2000.

F. Perrin, J. Pernier, O. Bertrand, and J. F. Echallier, “Spherical splines for scalp potential and current density mapping,” Electroencephalography and Clinical Neurophysiology, vol. 72, no. 2, pp. 184–187, 1989.

M. Hamedi, S.-H. Salleh, and A. M. Noor, “Electroencephalographic motor imagery brain connectivity analysis for BCI: a review,” Neural Computation, vol. 28, no. 6, pp. 999–1041, 2016.

J. Ginter Jr., K. J. Blinowska, M. Kamiński, P. J. Durka, G. Pfurtscheller, and C. Neuper, “Propagation of EEG activity in the beta and gamma band during movement imagery in humans,” Methods of Information in Medicine, vol. 44, no. 1, pp. 106–113, 2005.

P. Avanzini, M. Fabbri-Destro, R. Dalla Volta, E. Daprati, G. Rizzolatti, and G. Cantalupo, “The dynamics of sensorimotor cortical oscillations during the observation of hand movements: an EEG study,” PLoS One, vol. 7, no. 5, article e37534, 2012.

K. Sakihara and M. Inagaki, “Mu rhythm desynchronization by tongue thrust observation,” Frontiers in Human Neuroscience, vol. 9, p. 501, 2015.

R. Grech, T. Cassar, J. Muscat et al., “Review on solving the inverse problem in EEG source analysis,” Journal of NeuroEngineering and Rehabilitation, vol. 5, no. 1, p. 25, 2008.

B. He, Y. Dai, L. Astolfi, F. Babiloni, H. Yuan, and L. Yang, “eConnectome: a MATLAB toolbox for mapping and imaging of brain functional connectivity,” Journal of Neuroscience Methods, vol. 195, no. 2, pp. 261–269, 2011.

F. Babiloni, F. Cincotti, C. Babiloni et al., “Estimation of the cortical functional connectivity with the multimodal integration of high-resolution EEG and fMRI data by directed transfer function,” NeuroImage, vol. 24, no. 1, pp. 118–131, 2005.

A. C. Evans, D. L. Collins, S. R. Mills, E. D. Brown, R. L. Kelly, and T. M. Peters, “3D statistical neuroanatomical models from 305 MRI volumes,” in 1993 IEEE Conference Record Nuclear Science Symposium and Medical Imaging Conference, pp. 1813–1817, San Francisco, CA, USA, 1993.

G. Buzsáki, C. A. Anastassiou, and C. Koch, “The origin of extracellular fields and currents—EEG, ECoG, LFP and spikes,” Nature Reviews Neuroscience, vol. 13, no. 6, pp. 407–420, 2012.

H. Becker, L. Albera, P. Comon, R. Gribonval, F. Wendling, and I. Merlet, “Brain-source imaging: from sparse to tensor models,” IEEE Signal Processing Magazine, vol. 32, no. 6, pp. 100–112, 2015.

P. C. Hansen, T. K. Jensen, and G. Rodriguez, “An adaptive pruning algorithm for the discrete L-curve criterion,” Journal of Computational and Applied Mathematics, vol. 198, no. 2, pp. 483–492, 2007.

P. C. Hansen, “Regularization tools version 4.0 for Matlab 7.3,” Numerical Algorithms, vol. 46, no. 2, pp. 189–194, 2007.

C. Wilke, L. Ding, and B. He, “Estimation of time-varying connectivity patterns through the use of an adaptive directed transfer function,” IEEE Transactions on Biomedical Engineering, vol. 55, no. 11, pp. 2557–2564, 2008.

D. Kugiumtzis, “On the reliability of the surrogate data test for nonlinearity in the analysis of noisy time series,” International Journal of Bifurcation and Chaos, vol. 11, no. 7, pp. 1881–1896, 2001.

D. Kugiumtzis, “Surrogate data test on time series,” in Modelling and Forecasting Financial Data, pp. 267–282, Springer, Boston, MA, USA, 2002.

M. Kamiński, M. Ding, W. A. Truccolo, and S. L. Bressler, “Evaluating causal relations in neural systems: granger causality, directed transfer function and statistical assessment of significance,” Biological Cybernetics, vol. 85, no. 2, pp. 145–157, 2001.

C. W. J. Granger, “Investigating causal relations by econometric models and cross-spectral methods,” Econometrica, vol. 37, no. 3, p. 424, 1969.

T. Schneider and A. Neumaier, “Algorithm 808: ARfit—a matlab package for the estimation of parameters and eigenmodes of multivariate autoregressive models,” ACM Transactions on Mathematical Software, vol. 27, no. 1, pp. 58–65, 2001.

A. Delorme, T. Mullen, C. Kothe et al., “EEGLAB, SIFT, NFT, BCILAB, and ERICA: new tools for advanced EEG processing,” Computational Intelligence and Neuroscience, vol. 2011, Article ID 130714, 12 pages, 2011.

B. H. Jansen, J. R. Bourne, and J. W. Ward, “Spectral decomposition of EEG intervals using Walsh and Fourier transforms,” IEEE Transactions on Biomedical Engineering, vol. 28, no. 12, pp. 836–838, 1981.

G. Florian and G. Pfurtscheller, “Dynamic spectral analysis of event-related EEG data,” Electroencephalography and Clinical Neurophysiology, vol. 95, no. 5, pp. 393–396, 1995.

M. Sabate, C. Llanos, E. Enriquez, B. Gonzalez, and M. Rodriguez, “Fast modulation of alpha activity during visual processing and motor control,” Neuroscience, vol. 189, pp. 236–249, 2011.

M. Rubinov and O. Sporns, “Complex network measures of brain connectivity: uses and interpretations,” NeuroImage, vol. 52, no. 3, pp. 1059–1069, 2010.

O. Sporns, “Graph theory methods for the analysis of neural connectivity patterns,” in Neuroscience Databases: a Practical Guide, R. Kötter, Ed., pp. 171–185, Springer, Boston, MA, USA, 2003.

O. Sporns and R. Kötter, “Motifs in brain networks,” PLoS Biology, vol. 2, no. 11, article e369, 2004.

R. W. Floyd, “Algorithm 97: shortest path,” Communications of the ACM, vol. 5, no. 6, p. 345, 1962.

J. Goni, M. P. van den Heuvel, A. Avena-Koenigsberger et al., “Resting-brain functional connectivity predicted by analytic measures of network communication,” Proceedings of the National Academy of Sciences, vol. 111, no. 2, pp. 833–838, 2014.

J.-P. Onnela, J. Saramäki, J. Kertész, and K. Kaski, “Intensity and coherence of motifs in weighted complex networks,” Physical Review E, vol. 71, no. 6, article 065103, 2005.

G. Fagiolo, “Clustering in complex directed networks,” Physical Review E, vol. 76, no. 2, article 026107, 2007.

D. J. Watts and S. H. Strogatz, “Collective dynamics of ‘small-world’ networks,” Nature, vol. 393, no. 6684, pp. 440–442, 1998.

S. Achard, R. Salvador, B. Whitcher, J. Suckling, and E. Bullmore, “A resilient, low-frequency, small-world human brain functional network with highly connected association cortical hubs,” The Journal of Neuroscience, vol. 26, no. 1, pp. 63–72, 2006.

M. D. Humphries and K. Gurney, “Network ‘small-world-ness’: a quantitative method for determining canonical network equivalence,” PLoS One, vol. 3, no. 4, article e0002051, 2008.

O. Sporns, C. J. Honey, and R. Kötter, “Identification and classification of hubs in brain networks,” PLoS One, vol. 2, no. 10, article e1049, 2007.

A. Zalesky, A. Fornito, and E. T. Bullmore, “Network-based statistic: identifying differences in brain networks,” NeuroImage, vol. 53, no. 4, pp. 1197–1207, 2010.

Y. Benjamini and Y. Hochberg, “Controlling the false discovery rate: a practical and powerful approach to multiple testing,” Journal of the Royal Statistical Society: Series B (Methodological), vol. 57, no. 1, pp. 289–300, 1995.

M. Xia, J. Wang, and Y. He, “BrainNet Viewer: a network visualization tool for human brain connectomics,” PLoS One, vol. 8, no. 7, article e68910, 2013.

D. Cramer, Fundamental Statistics for Social Research: Step-By-Step Calculations and Computer Techniques Using SPSS for Windows, Routledge, New York, NY, USA, 1998.

D. Cramer and D. Howitt, The SAGE Dictionary of Statistics, SAGE Publications Ltd, London, UK, 2004.

D. P. Doane and L. E. Seward, “Measuring skewness : a forgotten statistic?” Journal of Statistics Education, vol. 19, no. 2, 2011.

N. M. Razali, Y. B. Wah, and M. Sciences, “Power comparisons of Shapiro-Wilk, Kolmogorov-Smirnov, Lilliefors and Anderson-Darling tests,” Journal of Statistical Modeling and Analytics, vol. 2, no. 1, pp. 21–33, 2011.

S. S. Shapiro and M. B. Wilk, “An analysis of variance test for normality (complete samples),” Biometrika, vol. 52, no. 3-4, pp. 591–611, 1965.

N. Kordjazi, A. Koravand, and H. Sveistrup, “Enhancing the representational similarity between execution and imagination of movement using network-based brain computer interfacing,” bioRxiv, 2017, In press.

N. Birbaumer and P. Sauseng, “Brain–computer interface in neurorehabilitation,” in Brain-Computer Interfaces, pp. 155–169, Springer, Berlin, Heidelberg, 2009.

F. De Vico Fallani, V. Latora, L. Astolfi et al., “Persistent patterns of interconnection in time-varying cortical networks estimated from high-resolution EEG recordings in humans during a simple motor act,” Journal of Physics A: Mathematical and Theoretical, vol. 41, no. 22, article 224014, 2008.

A. Athanasiou, C. Lithari, K. Kalogianni, M. A. Klados, and P. D. Bamidis, “Source detection and functional connectivity of the sensorimotor cortex during actual and imaginary limb movement: a preliminary study on the implementation of econnectome in motor imagery protocols,” Advances in Human-Computer Interaction, vol. 2012, Article ID 127627, 10 pages, 2012.

B. S. Baxter, B. J. Edelman, A. Sohrabpour, and B. He, “Anodal transcranial direct current stimulation increases bilateral directed brain connectivity during motor-imagery based brain-computer interface control,” Frontiers in Neuroscience, vol. 11, p. 691, 2017.

M. Grosse-wentrup, “Understanding brain connectivity patterns during motor imagery for brain-computer interfacing,” in Advances in Neural Information Processing Systems 21 (NIPS 2008), D. Koller, D. Schuurmans, Y. Bengio, and L. Bottou, Eds., pp. 561–568, Curran Associates, Inc., Red Hook, NY, USA, 2009.

Y. K. Kim, E. Park, A. Lee, C.-H. Im, and Y.-H. Kim, “Changes in network connectivity during motor imagery and execution,” PLoS One, vol. 13, no. 1, article e0190715, 2018.

W. Yi, S. Qiu, K. Wang et al., “Evaluation of EEG oscillatory patterns and cognitive process during simple and compound limb motor imagery,” PLoS One, vol. 9, no. 12, article e114853, 2014.

M. Hamedi, S.-H. Salleh, S. B. Samdin, and A. M. Noor, “Motor imagery brain functional connectivity analysis via coherence,” in 2015 IEEE International Conference on Signal and Image Processing Applications (ICSIPA), pp. 269–273, Kuala Lumpur, Malaysia, 2015.

S. B. Samdin, C.-M. Ting, S.-H. Salleh, M. Hamedi, and A. B. M. Noor, “Estimating dynamic cortical connectivity from motor imagery EEG using KALMAN smoother & EM algorithm,” in 2014 IEEE Workshop on Statistical Signal Processing (SSP), pp. 181–184, Gold Coast, VIC, Australia, 2014.

B. K. Kang, J. S. Kim, S. Ryun, and C. K. Chung, “Prediction of movement intention using connectivity within motor-related network: an electrocorticography study,” PLoS One, vol. 13, no. 1, article e0191480, 2018.

A. Athanasiou, N. Foroglou, K. Polyzoidis, and P. D. Bamidis, “Graph analysis of sensorimotor cortex functional networks – comparison of alpha vs beta rhythm in motor imagery and execution,” in SAN/NIHC 2014 Meeting, pp. 64-65, Utrecht, Netherlands, 2014.

R. B. Willemse, J. C. de Munck, J. P. A. Verbunt et al., “Topographical organization of mu and beta band activity associated with hand and foot movements in patients with perirolandic lesions,” The Open Neuroimaging Journal, vol. 4, pp. 93–99, 2010.

F. De Vico Fallani, F. Pichiorri, G. Morone et al., “Multiscale topological properties of functional brain networks during motor imagery after stroke,” NeuroImage, vol. 83, pp. 438–449, 2013.

A. M. Bastos and J.-M. Schoffelen, “A tutorial review of functional connectivity analysis methods and their interpretational pitfalls,” Frontiers in Systems Neuroscience, vol. 9, p. 175, 2016.

B. He, Ed., Modeling and Imaging of Bioelectrical Activity, Springer, Boston, MA, USA, 2005.

P. L. Nunez and R. Srinivasan, Electric Fields of the Brain: the Neurophysics of EEG, Oxford University Press, New York, NY, USA, 2006.