Reorganization of Resting-State EEG Functional Connectivity Patterns in Children with Cerebral Palsy Following a Motor Imagery Virtual-Reality Intervention

Motor imagery (MI) has been suggested to provide additional benefits when included in traditional approaches of physical therapy for children with cerebral palsy (CP). Regardless, little is understood about the underlying neurological substrates that might justify its supposed benefits. In this work, we studied resting-state (RS) electroencephalography (EEG) recordings of five children with CP that underwent a MI virtual-reality (VR) intervention. Our aim was to explore functional connectivity (FC) patterns alterations following this intervention through the formalism of graph theory, performing both group and subject-specific analyses. We found that FC patterns were more consistent across subjects prior to the MI-VR intervention, shifting along the anterior-posterior axis, post-intervention, for the β and γ bands. Additionally, group FC patterns were not found for the α range. Furthermore, intra-subject analyses reinforced the existence of large inter-subject variability and the need for a careful exploration of individual pattern alterations. Such patterns also hinted at a dependency between short-term functional plasticity mechanisms and the EEG frequency bands. Although our sample size is small, we provide a longitudinal analysis framework that can be replicated in future studies, especially at the group level, and whose foundation can be easily extended to verify the validity of our hypotheses.


Introduction
Cerebral palsy (CP) is the most common motor disability in childhood [1]. The term describes "a group of permanent disorders of the development of movement and posture, causing activity limitation, that are attributed to nonprogressive disturbances that occurred in the developing fetal or infant brain. The motor disorders of cerebral palsy are often accompanied by disturbances of sensation, perception, cognition, communication, and behavior; by epilepsy, and by secondary musculoskeletal problems" [2]. Since childhood is generally a highly physically active phase, CP can substantially impact a child's physical and cognitive development, affecting the quality of life of the children themselves, and of their primary caregivers and families [3][4][5].
Several treatments have been proposed over the years to maximize overall quality of life of CP patients; however, many of them display conflicting evidence: in 2013, Novak and colleagues reviewed 166 articles and identified 64 types of discrete interventions and interface and were able to achieve the experiment's goals. Despite the small number of children in this study (three children with CP), its results provided strong evidences for the feasibility of designing MI-based systems for children with CP. Results also reinforced the concept that each MI-related mental patterns are subject-specific and should be considered as such.
More recently, another study assessed the effects of associating MI with traditional physical therapy on upper limbs rehabilitation in children with CP. The authors measured patients' performance through the Assessing Hand Assessment (AHA) pre-and postintervention, and at a follow-up session. Participants were divided into a control group, which only received traditional physical therapy, and an experimental group, that received both the traditional and MI interventions. Overall, the authors reported that improvements in the AHA scores were larger for the experimental group, which were also mainted at the follow-up period [31].
Although the aforementioned results are certainly promising, no investigation regarding functional and/or structural neurological changes was carried out. Improvements are only measured through clinical scales and, despite being a valid and indeed useful tool for evaluating MI subjective performance and motor execution, they are limited in the sense that no quantitative information concerning the brain functioning throughout the therapy can be extracted. It could be expected that observed improvements would yield specific neural alterations that could be assessed with neuroimaging techniques. As a matter of a fact, several functional magnetic resonance imaging (fMRI) studies seeking neurological correlates of MI learning and practice often investigate alterations in resting-state (RS) (e.g., [32][33][34][35]), as this state provides the understanding of cortical interactions of the brain at a basal level (i.e., spontaneous brain activity) and could supply indicatives of neural plasticity [36]. When further including functional connectivity (FC), this type of analysis provides information regarding the brain's regional interaction mechanisms at rest, as well as their time evolution throughout the interventions under investigation. Despite being a prominent and promising analysis approach, to the best of our knowledge, studies of this type in CP patients are scarce.
With this perspective, in this work our aim was to propose an initial analysis framework for exploring FC patterns in resting-state EEG data of CP patients. We applied our methodology to a database of five children with CP in a similar VR intervention as described in [30]. We illustrate how our procedure can be used to describe alterations in said FC patterns comparing the resting-state activity before, and following, a combined action observation (AO) and MI VR intervention. To achieve this, we modeled the brain's functional interactions through graph theory, since employing this formalism establishes a mathematical representation of functional interactions and allows us to study topological properties of brain functional networks quantitatively. More specifically, we characterized the reorganization of functional links pre-and post-intervention at the group level, and described alterations in graph metrics for each patient. Our proposed data analysis procedure can be easily replicated, and even expanded upon, in future studies of the field.
The interface consisted in a virtual environment designed in Unity 2019.3.0a6 (Unity Technologies, San Francisco, CA, USA) ( Figure 1). Subjects followed a first-person avatar walking through a pre-defined path with a series of 22 obstacles, yielding a total session time of 18 min. When the avatar faced an obstacle, the walking stopped, and the patients should then relax, establishing a baseline period of about 4 s (times randomly varied between 3 s and 5 s in each obstacle to avoid stimulus habituation). Then, the obstacle was removed from their path, and participants were instructed to imagine the walking restarting, to resume it (motor imagery period; randomly varied between 3 s and 4 s). Patients were not actually controlling the interface, although they believed their neural signals were, indeed, commanding the virtual environment, whereas in fact the avatar's walking was restarted automatically.
Two RS-EEG recordings were made immediately prior, and post, to the AO+MI-VR intervention, and are the focus of the data analyses of this study. As briefly mentioned in the Introduction, this should enable the investigation of FC patterns related to spontaneous baseline brain activity. These recordings' duration was not uniform, with the average recording time being of (72 ± 15) s. Regardless, as explained below (Section 2.3), our approach takes relative thresholds for the data of each subject, which should minimize biases due to recording lengths.

Obstacle encounter
Obstacle removal Inter-blocks period Patients' demographic data are shown in Table 1, with information concerning their topographical classification-diplegia or tetraplegia-whether they had a premature birth, participants' age and their score in the Gross Motor Functional Classification System (GMFCS) [37].
Following that, a common average reference (CAR) filter, referencing electrodes to their average signal, was used with two goals. First, it should aid in reducing common artifacts present at all channels at the same time [38]. Additionally, CAR filtering has been suggested to, at least partially, aid in attenuating volume conduction effects [39]-since it theoretically removes a "common background signal" for all electrodes' time series-and is free of model-bound assumptions.

Graph Modelling of Functional Brain Networks
The methodological steps regarding the graph modelling of the brain are summarized in Figure 2, and are referred to throughout their explanation.

CORRENTROPY MATRICES FOR EACH SUBJECT AND TIME WINDOW -Vw[0]
... Illustrative diagram of the methodological steps: the recorded EEG time-series are analyzed in one-second windows, with the corresponding correntropy between sensors being computed. Next, subject-specific adjacency matrices are built based on binarization thresholds related to the strength of functional links and to the their persistence across the data recordings. Finally, our analyses encompass two approaches: one, at the group level, aims at determining changes in FC patterns common to all subjects; and the other, at the participant level, has as goal to characterize changes in the graphs' topology through specific metrics.

Graphs Adjacency Matrices
FC patterns were estimated using the correntropy with a Gaussian kernel [40]. The correntropy can be regarded as a generalized correlation function that can account for nonlinearities in the recorded time series [40]. Hence, since the brain can be considered as a non-linear dynamic system [41], and since EEG exhibits important non-linear characteristics [42], we chose to use this measure for computing similarities between the recorded time series of each pair of electrodes. Moreover, a recent application of correntropy to study classification scenarios of hands MI tasks showed that this measure was able to capture the expected spatial patterns of MI, while also performing better in terms of classification, when compared to other types of correlation metrics (namely, Pearson and Spearman correlations) [43]. Furthermore, even though other nonlinear approaches have been proposed to deal with these issues in the EEG signal, correntropy can be directly applied to the recorded data and, thus, generally becomes simpler and easier to compute [44]. We briefly describe the steps in the correntropy calculation in our study, below. A more through and mathematically formal discussion on correntropy can be found in [40,44]-the latter especifically concerning its application to EEG FC studies.
Consider the recorded EEG signals represented by , with a total of d channels and N data points, for each sample t in time window w. Let us omit the indice w in the quantity X (w) (t), for notation simplicity, and define a correntropy matrix V w [τ] for a given lag τ as such [40,43]: In Equation (1), the operator E[κ(x, y)] denotes the expected value over the kernel i,j (τ) denotes the correntropy coefficient between sensors i and j, from a total of d sensors, for a specific lag value. As mentioned, we worked with a Gaussian kernel [40], as follows: with σ being an adjustable parameter related to the kernel size. The kernel function produces a mapping of each vector X(n) − X(n − τ) over a sphere in the feature space and, in practice, its application computes a distance between the two vectors [40]. In our study, This would mean that maximum similarity of 1 is found when considering time series of the same electrode.
Values for V w [τ] in Equation (1) can be estimated through [40,43]: which is simply the mean of the application of the kernel function over all samples n, within the studied time window, for lag τ. To remain in the domain of functional connectivity, as no information regarding directionality of interactions and causality was needed (which would implicate effective connectivity instead [45]), we set τ = 0. Then, the crosscorrentropy matrix V w [0] represents the similarity matrix for data points contained in each time window w, between each pair of electrodes i and j: Estimations of V w [0] were performed for five frequency bands-α 1 , α 2 , β 1 , β 2 and γand over non-overlapping one-second time windows (that is, each time window can be represented by a matrix X(t) previously defined, with N = 250 data points). In the following, graph adjacency matrices were computed for each subject S-A S (ρ, r) = {a S (ρ, r) i,j } 32×32and were further binarized according to two thresholds: (i) ρ, which indicates the fraction of the top-strongest functional connections that should be maintained (i.e., not regarded as noise, or spurious); and (ii) r, which indicates the minimum fraction of time (from the total number of time windows) a link should be present to be considered. Thus, for example, setting ρ = 0.30 and r = 0.75 implies that only the top 30%-strongest connections that were present for, at least, 75% of the recording time, would compose that subject's adjacency matrix, A S (0.30, 0.75).
Nonetheless, the choice of thresholds for defining the final adjacency matrix is still an open problem in the literature [46]. Hence, for this study, we empirically defined an initial approach with a small set of threshold values; more specifically, we constrained: ρ = {0.3, 0.4, 0.5} and r = {0.55, 0.65, 0.75}. Such sets were considered neither too restrictive, which would largely limit the available number of functional links, nor too permissive, which could allow a high number of spurious FC links. Finally, employing relative thresholds ensures that all subjects will display the same number of functional connections in their adjacency matrix, and that differences in acquisition times will not compromise group analyses.

Types of Analysis
Analyses were performed at two levels, always comparing the pre-and post-session RS-EEG recordings: (a) a group level, considering an adjacency matrix describing only the common FC links for all subjects; and (b) the individual level, contemplating the distribution of graph metrics values for each subject.
The group level analysis considers a group graph-G(ρ, r) = {g ij (ρ, r)} 32x32 -accounting for only the functional connections present for all subjects for a given (ρ, r) pair. Thus, given M subjects, each element of G(ρ, r) is defined as: Given the relatively small total number of links for G(ρ, r), only analyses regarding the links themselves were performed (i.e., we considered the group adjacency matrices to be too sparse to yield meaningful graph metrics).
The intra-subject analysis involved exploring graph metrics' distribution over each electrode during the pre-and post-intervention conditions (since the larger number of links, in turn, complicates a similar analysis approach to the group study). To achieve this, we compared the metrics' distributions, over time, during each condition (i.e., preand post-session) and tested for statistically significant differences through a Wilcoxon Sign-Rank test (correcting the p-value significance threshold of 0.05 with the Bonferroni correction, considering 5 subjects, 4 graph metrics, 32 electrodes and 5 frequency bands). We then explored which regions displayed statistically significant increasing and decreasing trends for each patient. Four metrics were analyzed: (1) the degree, (2) the clustering coefficient, (3) the node path length, and (4) the betweenness centrality. All measures were estimated through the routines of the Graph Connectivity Toolbox for MATLAB, by Rubinov and Sporns [47]. Below we briefly describe the mathematical meaning of these metrics, and their neurophysiological interpretation in our study is left to be discussed in the appropriate section.
A graph's degree at node i (k i ) simply expresses the number of connections it makes with the rest of the network, thus being mathematically expressed as the sum of the binarized adjacency matrix's elements: A node's clustering coefficient (CC i ) considers mainly the connections among its closest neighbors, and can be seen as a measure of how many closed triplets can be formed, in relation to the possible total that could exist, around that node. Thus, a node with high CC has its nearest-neighbors making a large number of connections among themselves, and vice versa. It can be expressed as: The path length between nodes i and j-l ij -, for a binary graph, measures the least amount of edges to be travelled between these two nodes. The average path length for node i, l i , is estimated through averaging l ij .
Finally, intuitively, the betweenness centrality (BC) measures how often a node acts as a "bridge" between any two other nodes in the graph. More formally, the BC for node i can be defined as the number of shortest path lengths to reach node l, starting from node j, that necessarily pass through node i-σ jl (i), in relation to the total number of shortest paths between j and l (that may not involve i)-σ jl : The factor 2 (N−1)(N−2) normalizes BC within [0,1]. Figure 3 displays the number of FC links for G(ρ, r) according to the studied frequency bands and thresholds. Overall, regardless of the thresholds' choice, the same tendency can be observed for the group's adjacency matrix, by means of the number of functional links shown in Figure 4, which tended to increase (decrease) for the higher (lower) frequency ranges. Indeed, results in Figure 3 indicate that FC patterns are more consistent across subjects for the higher frequency bands, with no common functional links being found in the α range. Furthermore, results in Figure 3 also indicate that, with exception of β 1 for some (ρ, r) combinations, FC patterns were more consistent across subjects prior to the intervention.   Figure 4 displays the FC links for G(ρ, r), for the three frequency bands that presented a non-empty group adjacency matrix, for both the pre-and post-session conditions. Following the color scheme in Figure 3, functional connections present only during the pre-intervention recording are shown in blue, whereas the ones found only at the postsession recording are displayed in red. Connections that were present at both instances are shown in green. These results allow complementing the findings in Figure 3 by providing a notion that the functional reorganization following the AO+MI-VR practice not only generally results in a lesser number of common links, but it also shifts their locations to more posterior nodes.

Individual Results-Study of the Networks' Topology through Graph Metrics
Despite the aforementioned results for general trends at the group level, each brain is unique and, especially in the case of neurologically damaged patients, brain plasticity and functional reorganization mechanisms can be especially important to guide intervention protocols [48]. Therefore, intra-subject analyses were also performed, to explore possible subject-specific effects of the AO+MI-VR intervention. The following results were obtained setting the thresholds values to the most restrictive possibility explored in this study; that is: (ρ, r) = (0.30, 0.75). For all metrics, variations were simply taken as the difference between their values at the corresponding nodes between the post-and pre-intervention periods. Hence, positive variations indicate increases in the metric following the AO+MI-VR session, and vice versa. We emphasize that these analyses are highly exploratory. Nevertheless, we attempted to extract common patterns, whenever possible.
Degree: Figure 5 shows results for the degree for each subject (rows), and for all frequency bands (columns). The circles' radii around each electrode site directly relates to the magnitude of the metric's change, as also represented by the colors scale.
The main sites for degree changes ( Figure 5) are predominantly on the areas indicated by the group's graph-i.e., frontal and parietal sites. Variations on the C-labeled electrodes (except for S03 in the α range and S01 in β 2 and γ) could also be observed, although to a lesser extent.
Clustering Coefficient: Figure 6 displays results for the CC. Empty plots should be interpreted as not containing statistically significant variations. Overall, similarly than for the degree (Figure 5), alterations in the CC were more often found over frontal and parietal nodes. However, in this case, most prominent variations were especially localized closer to the scalp's midline-which is especially true for S01. This means that, even though the total number of connections generally varies all over the scalp ( Figure 5), alterations in the CC are, at least to some extent, more restricted to sites closer to the scalp's midline. Regardless, as with the degree, the type of change (i.e., increase or decrease) is highly subject-dependent.
Node path length: Even though alterations in the degree and clustering coefficient were observed, no statistical significant changes were found for the node path length with our methodology.
Betweenness Centrality: Figure 7 shows results for the BC. Differently than for the other metrics, and with only a few exceptions (namely, S05 in the α range), the BC generally displayed, most noticeably, decreases. In addition, nodes with increased BC overall present a smaller magnitude (which can be seen through the smaller circle's radii). Once more, results are largely subject-specific and emphasize frontal and, parietal nodes, with less stress on central electrodes.

Group Analysis
Given the definition of the group graph, G(ρ, r), a larger number of links implicates that the FC patterns were found more consistently across subjects, and vice versa. Evidently, less restrictive thresholds should yield larger numbers of functional connections; however, thresholds that are too low have a greater risk of resulting in spurious links and, hence, allowing for misleading information. Nonetheless, as mentioned, there is still no established manner for setting thresholds for this type of study. Therefore, by working with different sets of values, we should be able to assess whether general trends can be observed at the group level, at least regarding specific thresholds ranges. It is worth emphasizing that, given the definition of the elements of G(ρ, r) in Equation (5), observed trends in Figure 3 necessarily occurred for all patients. Since distinct EEG rhythms might relate to different cognitive and affective states [49], it is actually reasonable that FC patterns were not equal across all rhythms. In addition, it is worth observing that, even though MI is predominantly linked to rhythms within the α range (i.e., the µ band), this task could also induce the observed alterations for the β and γ bands; the latter, especially, is commonly overlooked in MI-EEG studies (notwithstanding, recent works have increasingly argued for the involvement of this band in this type of task-e.g., [50]). Moreover, as mentioned, functional links were more consistently found across subjects prior to the AO+MI-VR intervention.
Interestingly, a functional near-infrared spectroscopy (fNIRS) study in children with CP [51] has showed that, immediately following constraint-induced movement therapy, FC patterns' consistency across subjects diminished between the supplementary motor area and the premotor and primary sensorimotor areas. However, the original pre-therapy consistency patterns were recovered at a 6-months follow-up period [51]. Although we have employed distinct analysis and MI intervention approaches, and the nature of the EEG and fNIRS signals certainly differ, our results in Figure 3 did display less consistent inter-subject FC patterns post-session, and, complementarily, Figure 4 provided the notion that the intra-subject consistency in functional links shifted from premotor and primary sensorimotor areas mainly towards the parietal cortex. In this paper, we obtained, at least to some extent and qualitatively, with a distinct methodology for data analysis, therapy intervention and brain activity recording technique, similar results to Cao and colleagues in [51].
Furthermore, a recent paper investigated the relationship between the frontoparietal network reconfiguration, assessed through EEG FC analysis, and attentional performance in healthy adults [52]. Overall, the authors propose that weak frontoparietal connectivity could be an indicator of improvement in performance of everyday tasks involving attentional performance due to the capacity for network reconfiguration involving these functional connections [52]. Therefore, the ability to learn was related to the degree of frontoparietal connections reorganization. In our study, indeed, frontal and parietal regions were the ones mostly implicated in the functional reorganization of the group's adjacency matrix in our study ( Figure 4).
Moreover, the frontoparietal network has been consistently implicated in higher cognitive control and considered as a functional hub that interacts with other brain networks, being flexible in its functional communications to allow the learning of new tasks involving cognitive control [53]. In relation to MI, magnetic resonance imaging (MRI) works have also indicated short-term plasticity mechanisms involving both structural and functional alterations, with changes in functional connectivity involving the defaultmode-network, and emphasizing that FC changes do spread beyond primary sensorimotor networks [54,55]. Given that those studies involved MRI, rather than EEG, anew, at this point, the shift in the functional links observed in Figure 4 can only be conjectured to relate to the aforementioned characteristics of the frontoparietal network reorganization mechanisms and learning. Certainly, this requires further investigation by, for example, attempting to quantify the functional reorganization degree and correlate it to the patient's MI performance.
Moreover, Figure 4 also suggests a right-hemisphere preference for the β 2 and γ bands, being slightly less evident for β 1 at a few (ρ, r) combinations. At this point, however, we are unable to assess whether this FC links lateralization might relate to any brain functional mechanisms, or even to patient lesion sites.
Finally, the functional connections present at both recordings were predominantly located over primary sensorimotor cortex-related electrodes, with a particular clusterization surrounding C4, especially for the β 2 and γ bands. On the one hand, since short-term functional plasticity largely depends on the performed task [56], this might relate to the MI task performed by the patients. On the other hand, the common FC links present at both recordings might not be strictly motor-related, since they also involve, at least for some instances in Figure 3, functional interactions with areas on the posterior cortex and, thus, could be an effect of more cognitive facets of MI, or even a result of the AO activity.
In summary, our group results indicate that: (1) RS-EEG-FC patterns are more consistent across subjects prior to the AO+VR-MI intervention, being mostly predominant on higher frequency rhythms (β 1 , β 2 and γ) and located over frontal and central nodes; (2) immediately after the intervention, such patterns are less consistent across subjects, for all studied frequency bands, and they shift to more parietal areas, with a predominance for right-hemisphere functional interactions, slightly more evident for the β 2 and γ bands; (3) FC links common to both pre-and post-session patterns were mostly bound to the electrodes generally located over the primary sensorimotor cortex, with a particular clusterization around C4 (and even other electrodes, such as C6, CP4 and CP6). It remains to be investigated whether the degree of group reorganization correlates with overall task performance and if the observed changes are maintained in the long-term.

Intra-Subject Analysis
Degree: According to our methodology, functional connections ultimately represent the extent of similarity between two electrodes' time series and, thus, an increase (decrease) in the degree of the i-th node can be interpreted as an increased (decreased) overall similarity between time series recorded at i and the other neighbor nodes. Hence, even though subjects display specific sets of spatial alterations in FC patterns (analyzed here through the graph metrics variations), at least regarding the degree (i.e., the number of connections at each node), frontal and parietal sites do seem to play a critical role in the short-term plasticity (as also indicated by the group results). Nonetheless, a common pattern of synchronicity alteration could not be found for all subjects: for example, subjects S04 and S05 displayed degree increases on frontal sites in the α range, whereas S02 showed pronounced decreases. This only emphasizes each patient's uniqueness.
In addition, for all subjects, the graphs' topology alterations patterns are more distinct among the α and β 2 bands than between α 1 and α 2 , or between β 2 and γ bands; i.e., alteration patterns are more consistent both for within the α range, as well as within β 2 and γ. Therefore, we hypothesize that the observed short-term plasticity mechanisms, besides being frequency band-specific (as already suggested at the group level by Figures 3 and 4), seem to behave as if the lower (α 1 and α 2 ) and upper (β 2 and γ) frequency bounds constituted distinct states, or phases, of the brain's functional reorganization, with results for β 1 providing a "transition" frequency band between the two. In fact, the change in the degree can occur distinctly at the same electrodes for different frequency bands; that is: nodes that present increases in one band can show a decrease in another, and vice versa-such as P3 for S04 and FC2 for S05, supporting this hypothesis.
Clustering Coefficient: The CC for a node i considers mainly the connections among its closest neighbors. Thus, the neighbors of a node with high CC are largely connected among themselves, and vice versa. Even though, in a sense, CC also accounts for the number of connections, this metric reflects more local interactions than the degree. The concomitant interpretation of both metrics can, then, provide complementary information. A few common trends found for the degree can be observed for the clustering coefficient as well (Figure 6), such as the frequency-dependent patterns and the highly subject-specificity of results. Notwithstanding, for the CC, none of the patients showed significant variations for the α 1 graph.
Moreover, it is important to emphasize that, even though decreases in the total number of connections were mainly associated with a correspondent loss in the local number of links (i.e., concomitant decrease of degree and CC), this is not necessarily true for all instances, such as with FC2 for S04 at γ. This then illustrates that specific nodes can display a decrease in their total number of connections, even if their clusterization (emphasized by the nearest-neighbour connections) increases.
Node path length: Not finding any type of statistically significant result for the path length might be due to an ill-posed problem by isolated nodes in the network, possibly caused by thresholding, which theoretically sets their path length to mathematical infinity. We are still reporting it here, though, to stress the current existing need of research to address optimal threshold choices for graphs representing brain functional networks.
Betweenness Centrality: The evident decreases ( Figure 7) indicate a "loss of importance" in that specific node regarding information flow in the network. Hence, our results suggest that, overall, regulation of information flow tended to be less centralized at specific regions. Still, an exception to this was clearly observed for S05 in the α band. Note, though, that this general increasing trends shift to decreases, more consistently with the other patients, when transitioning from the α to the higher frequency ranges, which also supports our hypothesis for distinct plasticity mechanisms that are frequency-specific.
To summarize, investigation alterations in the graphs' topology at the individual level evidenced that: (1) as expected, alteration patterns are highly subject-specific. Even though common trends can be found at the group level, the group analysis by itself does not provide the whole picture for understanding the brain dynamics of each patient; (2) short-term plasticity mechanisms might be intrinsic to each frequency band.

Limitations and General Considerations
Variability in brain data can be extremely large and can pose several problems when seeking general or average patterns across both healthy subjects and neurological patients. However, it can also be very informative as it emphasizes that each brain is unique and subjects may present their own particularities in responding to an intervention, even when common patterns can indeed be found for a group of participants; in other words, brain variability can also provide relevant data, and should not be solely considered as noise [57]. With this under consideration, our analyses explored alteration patterns occurring for all subjects, as well as individually.
Since we are comparing the moments immediately prior, and post, the AO+MI-VR intervention, and amongst the outlined changes we also identify variations that occurred for all subjects in motor related areas, compatible with MI-related changes found in other works [32,58], our study assumes that detectable changes in RS-FC patterns are at least partially due to the AO+MI VR task that the patients performed. Additionally, frontoparietal or occipital changes might be due to the attention requirements or visual stimuli implied in task. However, the determination of the specific contributions would need additional research.
Based on our results, we are currently unable to attribute any differences between subjects to the available demographic factors, as no sub-grouping was evident with the accessible quantities. Nonetheless, general hypotheses that could be drawn from this initial study include the larger consistency in FC patterns, at rest, across subjects prior to the MI-VR intervention and the distinct short-term modifications mechanisms at play in these patterns according to different EEG frequency bands.
In any case, even these more general hypotheses still require further investigation in larger datasets, as our small number of participants and, as a result, the lack of statistical power, impairs the extrapolation of our conclusions. In addition, the neurological damage is patient-specific, which can further hinder generalized observations. Regardless, our methodology can be easily replicated in future studies and contributes in providing quantitative basis for the investigation of FC patterns in CP patients submitted to MI-based interventions, a research topic still lacking on this sort of analyses procedures.
Finally, despite our FC analysis being performed at the sensor level, instead of the source domain, a 2018 study showed that, at least regarding the global connectivity (i.e., the overall indices in the adjacency matrices), estimates of FC in both domains (sensor and source) were closely correlated [59]. They generally concluded that FC can be reliably estimated, too, from the sensor space. However, they do stress that analyses should be carefully considered, and that metrics that account for volume conduction effects when estimating FC should be preferred [59]. Regarding these points, the CAR filtering in our preprocessing steps should, at least partially, minimize volume conduction effects. Moreover, the correntropy, the measure we employed in our study for estimating the FC patterns, has been suggested to be robust to noise and able to detect non-linear relationships between EEG signals [44]. It has also suggested to be more robust than other more traditionally employed metrics for EEG-FC studies [44]. Besides, studies from our group have also evidenced correntropy as a solid alternative to other well-established correlation metrics, being able to provide additional information the others might miss [43].

Conclusions
We have provided an analysis framework for studying FC patterns in children with CP by presenting alterations at both the group level -indicating observed common trends for all patients -and at an individual level, which should provide more subject-particular insights. There is a considerable number of variables to be investigated, as well as several methodological choices that remain open problems in the literature.
Given the low number of common functional links for the group graph, only a study regarding the (re)organization of functional connections was conducted. In this case, our main findings were: (1) FC links are more consistent across subjects prior to the VR-MI task; (2) the reorganization between the pre-and post-session recordings shifts connections along the anterior-posterior axis; and (3) there is a lateralization of functional links on the right hemisphere. We raised the hypothesis that (1) and (2) might relate to MI-learning and short-term functional plasticity mechanisms. However, we currently have no speculations for (3).
Regarding the intra-subject analysis, as expected, alterations in the graphs' topology are largely subject-specific, even though a common "back-bone" of network alterations could be observed through the group's graph. Although results with this type of analysis permit a large array of interpretations that can be drawn for each subject, overall patterns showed dependence between graph metrics' variations and the frequency bands, suggesting that short-term functional reorganization in the brain could be frequency-specific. In particular, our findings hint at three frequency band-dependent facets of such reorganization mechanisms occurring for the lower (in our case, α 1 and α 2 ) and higher (β 2 and γ) frequency ranges, with an intermediate phase represented by β 1 .
Future analyses may expand on ours for a larger group of subjects to verify whether the observed general tendencies are maintained. Additionally, our analysis framework can be extended to other thresholds ranges and methodologies of FC estimation. This can lead to valuable and constructive comparisons that should aid in enhancing our current understanding on RS-EEG-FC patterns of children with CP, as well as on the influence of MI interventions on such patterns.  Data Availability Statement: Data will be shared upon reasonable request subjected to parental consent.