Brain electrical microstate features as biomarkers of a stable motor output

Objective. The aim of the present study was to elucidate the brain dynamics underlying the maintenance of a constant force level exerted during a visually guided isometric contraction task by optimizing a predictive multivariate model based on global and spectral brain dynamics features. Approach. Electroencephalography (EEG) was acquired in 18 subjects who were asked to press a bulb and maintain a constant force level, indicated by a bar on a screen. For intervals of 500 ms, we calculated an index of force stability as well as indices of brain dynamics: microstate metrics (duration, occurrence, global explained variance, directional predominance) and EEG spectral amplitudes in the theta, low alpha, high alpha and beta bands. We optimized a multivariate regression model (partial least square (PLS)) where the microstate features and the spectral amplitudes were the input variables and the indexes of force stability were the output variables. The issues related to the collinearity among the input variables and to the generalizability of the model were addressed using PLS in a nested cross-validation approach. Main results. The optimized PLS regression model reached a good generalizability and succeeded to show the predictive value of microstates and spectral features in inferring the stability of the exerted force. Longer duration and higher occurrence of microstates, associated with visual and executive control networks, corresponded to better contraction performances, in agreement with the role played by the visual system and executive control network for visuo-motor integration. Significance. A combination of microstate metrics and brain rhythm amplitudes could be considered as biomarkers of a stable visually guided motor output not only at a group level, but also at an individual level. Our results may play an important role for a better understanding of the motor control in single trials or in real-time applications as well as in the study of motor control.


Introduction
During action execution, the brain operates an adaptive control of movements by constantly comparing the incoming sensory signals with the predicted sensory outcome in order to successfully tune the motor output to the desired performance. During this process, the values of the parameters of the control rule on which motor control is based undergo changes due to the external sensory feedback received during motor action. Among the various external sensory feedbacks that can be provided to the motor system, the key source of information to achieve a proficient motor control is visual information [1,2]. This process of sensory feedback, integration and motor execution is supported by recurrent loops between parietal and motor areas within distributed fronto-parietal networks [3][4][5][6][7]. Accordingly, when maintaining a constant output, as in the case of the application of a constant force during an isometric muscle contraction, sensory feedback and integration processes form the foundation for a good performance, i.e. the maintenance of a required constant contraction level [8,9].
Although it is known that the cortical oscillatory activity, as measured by electroencephalography (EEG) or magnetoencephalography, plays a fundamental role in different mechanisms regulating both movement execution and cognitive processes [10], its functional role during the control of a stable motor output is not yet fully understood. Beta band (15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25) activity in the motor network has been traditionally associated with a steady-state motor output [11][12][13][14][15]. Indeed, a higher power in the beta band was found to be correlated to maintaining a more stable force level [12], and the decrease of the beta band amplitude (namely event related de-synchronization (ERD)) was typically found during the preparation, execution, imagination, or observation of movements [16]. Nevertheless, no correlation has been found between the modulation of the beta ERD during motor task execution and the values of several parameters characterizing the maintenance of a constant force, such as movement speed, weight of the manipulated load and level of the exerted force during an isometric contraction [17,18]. In 2010, Engel and Fries hypothesized that the beta band activity can represent the signature of an active process needed to maintain the actual motor or cognitive state by impeding the elaboration of new movements or cognitive processes [10]. Within this framework, evidence showed that the beta band is involved in a more efficient elaboration of the sensory feedbacks needed to monitor the state of the motor system [19]. To support this hypothesis, other studies found that the beta band activity can change in relation to the expectancy of an incoming event [20] while inhibiting the processing of new movements [21]. For these reasons, in the context of an adaptive control of movement, the beta band activity could play a controversial role in the maintenance of a constant motor output because of a sort of 'inertia' of the motor control system. In fact, the activity of the beta band appears to be related to the inhibition of unexpected external feedback in favor of the top-down endogenous control [10].
Beyond the functional role of cortical oscillations during a stable motor output, the brain dynamics associated with a successfully maintained isometric force is still unexplored. So far, studies investigating this issue [11,12,14] considered the brain electrical field only in specific brain areas with an a priori selection of electrodes of interest, at given time intervals and/or in specific frequency bands. However, to better understand the sensory feedback and integration processes during the control of a good isometric motor performance, methods estimating the global dynamics of widespread neural networks that do not rely on any type of a priori hypothesis (such as the electrode location or specific frequencies) are needed. Microstate analysis permits to represent the multichannel EEG signals in terms of global neuronal activity and to identify periods in which the topographic configuration of the electrical potential on the scalp is quasi-stable [22,23]. This approach describes the neuronal activity by a sequence of rapid transitions among a limited number of dominant scalp potential topographies or 'maps' , each lasting about 40-120 ms, which account up to about the 70%-80% of the total topographic variance [22]. The epochs of topographic quasi-stability, referred to as 'microstates' , can capture the activity of distributed cortical networks as global patterns of scalp potential topographies that dynamically vary over time in an organized manner [22,24]. Most previous studies converged on the existence of four dominant microstates in healthy adults at rest [22]. These four microstates exhibited highly similar topographies across studies and seemed to consistently dominate the EEG data across different age ranges, conditions, and pathological states [22,25,26].
The aim of the present study was to elucidate the brain dynamics during the maintenance of a steady force. We explored the correlation between the temporal evolution of subject-specific microstate features and the motor performance, defined in terms of maintenance of a fixed level of the exerted force in a visually guided isometric contraction task. Traditionally, microstates are extracted from broadband EEG signals and the link between the broadband topographies and brain activity in specific frequency bands remains unclear. Several studies have suggested that broadband topographies are not dominated by a specific narrow-band brain activity and that the microstate dynamics cannot be derived from one frequency band or a subset of specific EEG frequencies [27,28]. From this perspective, the quantification of microstate features and the modulation of brain rhythms are considered complementary as they could explain different aspects of motor control. For this reason, we also calculated the amplitude of brain oscillations in given frequency bands and their role during the maintenance of a constant force.
The search for a link between microstate features and brain rhythms during the maintenance of a constant contraction can be considered a multivariate regression problem [29]. In this type of problem, a high number of input features are included in the model (microstate features and brain rhythm amplitudes) to predict the output (motor performance). Our approach relied on a data-driven method [29,30] with the avoidance of restricting a priori assumptions. Specifically, because of a possible strong correlation between the input variables, a partial least square (PLS) regression analysis [29,31], that accounts for such correlations, was employed to reduce overfitting of the data and to produce more generalizable results.

Methods
The section 2 is structured as follows: we first describe the type of volunteers participating in the experimental studies (section 2.1), the type of experimental task employed (section 2.2) and how the motor performance of the volunteers was assessed and evaluated (section 2.3). Then we describe how the EEG data was recorded and analyzed (section 2.4). This subsection includes a specific section dedicated to describing the microstate analysis at rest and during isometric contraction (section 2.4.1) and the regression model used to verify the predictability of the motor performance using the EEG band amplitudes and the microstate metrics (section 2.4.2). Finally, there is a subsection dedicated to describing the statistical analysis performed (section 2.5).

Subjects
Eighteen healthy volunteers (twelve males, mean age of 21.3 ± 1.8 years) were recruited for the study. All subjects reported no history of neurological or psychiatric disease and did not receive psychoactive medications at the time of the recording. The Edinburgh questionnaire was used to ensure the righthandedness of all subjects [32].
The research was conducted in accordance with the principles of the Declaration of Helsinki. The experimental protocol was approved by the local Ethic Committee (ID number: richlhf92) and all subjects signed the informed consent form.

Experimental task
Subjects were requested to press a bulb positioned on the palm of their right hand by moving the right thumb against the other fingers, and to maintain a stable level of contraction. A pressure sensor was integrated in the bulb to measure the exerted force, which was represented on a screen by the position of a horizontal black bar that can slide on a vertical green bar (figure 1). The requested constant force level was indicated by two black lines at the sides of the vertical green bar. The maximal voluntary contraction (MVC), measured in terms of pressure with the pressure sensor, was acquired at the beginning of the experiment and expressed as the mean of three repetitions of maximal contractions, each one maintained for 1 s. The experimental task consisted of 48 trials. Each trial was composed of alternating resting state and contraction phases. The experiment started with a resting state of 6 s duration (figure 1). At the end of this period the contraction phase, lasting for a total of 8 s, started: two black lines appeared at the middle height of the vertical green bar-corresponding to a reference level of contraction equal to 15% of the MVC-to prompt the subjects to press the bulb, reach the indicated contraction level, and maintain it constant for 4 s. At the end of these 4 s, the two black lines moved to a different position, indicating a new level of contraction that the subject had to reach and maintain for 4 more seconds. At the end of this second contraction, the two black lines disappeared indicating that a new resting state phase of 6 s started (figure 1). During the second part of the contraction phase, the two black lines were set at one of four levels of contraction, two of which were below the reference level of contraction and corresponded to 5% or 10% of the MVC, whereas the other two were above the reference level of contraction and corresponded to 20% or 25% of the MVC. The four different levels of contraction were balanced in the experimental task (12 occurrences for each contraction level) and presented randomly to avoid learning effects.

Motor performance evaluation
The relative standard deviation of the exerted force with respect to the expected force level was calculated as an index of the stability of the exerted force during a given time interval. Specifically, if x 1 , x 2 , … x N are the N samples of the exerted force during the ith time interval (for the definition of the time intervals see section 2.4.1 below) and RFL i is the required force level for the same time interval, both expressed as percentage of MVC, the force variation index for the ith time interval (FV i ) was calculated as: The more stable the exerted force is during the contraction, the lower is the FV i value. An FV i value equal to zero corresponds to a level of the exerted force equal to the RFL i over the time interval, indicating a maximal stability, hence the best performance. On the contrary, a high FV i value corresponds to values of exerted force different from the RFL i , indicating a low stability and a poor performance.

EEG recording and data analysis
Brain electrical activity was recorded during the performance of the experimental task by means of a 128 channel EEG system (Net 300, Electrical Geodesics Inc., USA). Skin-electrode impedance was measured before each EEG recording and kept below 50 kΩ, which is the recommended impedance value below which EEG signals recorded with an Electrical Geodesics system can be considered of good quality (www.egi.com/knowledge-center/item/67the-effects-of-high-input-impedance-on-eeg-dataquality). EEG data were measured with a sampling frequency of 1000 Hz and processed offline. The exerted force was measured by the pressure applied on the pressure sensor with a sampling frequency of 100 Hz. The pressure sensor was connected to a digital board hosting in-house-developed electronics that recorded the exerted pressure while providing the visual feedback to the subjects and sending triggers to the EEG system for synchronization purposes (Interactive Pressure Sensor, InPresS [33]).
The raw EEG data were filtered with a 2nd order forward-backward Butterworth bandpass filter with cut-off frequencies at 0.5 Hz and 45 Hz. Filtered EEG data were visually inspected to exclude EEG channels exhibiting excessive noise, poor scalp-surface contact, or isoelectric saturation [34]. Excluded EEG channels were replaced by spline interpolation of the neighboring channels [35]. The EEG trials with more than 50% of the electrodes exhibiting excessive noise were trimmed from the data. A semiautomatic procedure, based on independent component analysis, was applied to remove ocular, cardiac or electromyogenic artifacts [34,35]: filtered EEG data were pre-whitened by principal components analysis and decomposed into 50 ICs using the extended Infomax algorithm, which can better separate signal sources that may exhibit super-Gaussian and sub-Gaussian distributions [36]. An average number of 9 ± 4 artefactual components (mean across subjects ± standard deviation) were identified for each subject. The artefactual ICs were disregarded, and the clean EEG signals were reconstructed by reprojecting the retained non-artefactual ICs onto the sensor space.

Microstate analysis at rest and during isometric contraction
To assess possible differences in microstate templates and dynamics between rest and movement intervals, the microstate templates were separately extracted during the 3 s time intervals of rest preceding the onset of the first contraction and during the last 3 s of the second contraction (from 5 to 8 s of the total contraction phase). The 3 s of the second contraction correspond to the force levels exerted at 5%, 10%, 20% or 25% MVC. We chose to analyze the last 3 s of the second contraction because this interval is well separated from the rest period and therefore it is not influenced by transient phenomena due to the rest/ movement adaptation. Moreover, the time interval between 5 and 8 s from the movement onset corresponds to a period of high stability, in which the subject was accurately maintaining the new requested force level (5%, 10%, 20% or 25% of MVC). The two subsequent force levels have been introduced in the experimental task to simulate conditions of contraction modulation as they occur in ecological motor performance. Finally, the different levels of the exerted force across trials (from 5% to 25%) simulate the variability in motor performance, needed for subsequent regression analysis.
To perform the microstate analysis, the EEG signals during the intervals defined above were concatenated separately for the rest and contraction conditions. EEG data were referred to the common average. The global field power (GFP), defined as the standard deviation across channels, was calculated. The spatial topographies corresponding to the maximum values of the GFP were then submitted to the k-means algorithm, a modified version of a clustering algorithm [37]. The optimal number of clusters was chosen by repeating the k-means algorithm with k varying from 1 to 12 and estimating the Krzanowski-Lai (KL) criterion for each repetition. The number of clusters corresponding to the second KL maximum value was chosen as the optimal number of clusters [24]. This procedure returns k templates for each subject for both rest and movement conditions. To obtain global templates across conditions the mean of the topographies (across subjects) were separately calculated at rest and during movement by applying a second clustering procedure. With this approach, k global topographic templates for rest and k global topographic templates for movement were obtained. The match of the global topographic templates between rest and movement conditions was based on the minimum of the dissimilarity value between the templates of the two conditions. The global dissimilarity between pairs of global templates is calculated as: where u i and v i are the electrical potentials of the ith electrode in the global topographic templates u and v respectively (the global topographic templates u and v refer to the same microstate template in the rest and movement conditions, respectively); GFP u _ and GFP v _ are the GFPs of the microstate templates u and v; N is the number of electrodes. Rest and movement topographies were compared using the topographical analysis of variance (TANOVA, [38]), separately for each microstate template pairs. Microstate sequences were then computed by fitting the templates corresponding to each condition (rest or movement) to the continuous original filtered EEG data. This back-fitting procedure was done by calculating the maximum spatial correlation between each global template and the topography at each time instant and assigning to each time instant the template with higher correlation (winner-takes-all procedure). Given that residual noise could affect the filtered EEG signals possibly producing unstable instantaneous topographies, a temporal smoothing procedure was applied ensuring the same global template is assigned to EEG intervals of at least 30 ms [39].
For each subject and for each global template, the following metrics were calculated for both rest and movement conditions, as mean values across trials: (a) Mean microstate duration: the mean of the time covered by a single microstate template, expressed in milliseconds. (b) Mean occurrence per second: mean number of distinct occurrences of a given global template within the time interval [40]. (c) Global explained variance: a metric expressing how well the global template describes the whole dataset [24].
We also calculated the transition probabilities for each pair of global templates: the X to Y transition was evaluated as the percentage of the number of transitions from microstate X to microstate Y over the total number of transitions in the considered time interval. Similarly, the Y to X transition was the number of transitions from microstate Y to microstate X over the total number of transitions in the considered time interval. Directional predominance quantifies the directional asymmetries in the transitions between two microstates and was calculated as the difference between the X to Y and Y to X transitions. All possible directional predominance values were calculated for all possible microstate pairs. A positive value of X↔Y indicates a higher number of transitions from X to Y than from Y to X, whereas a negative value indicates the opposite predominance of transitions from Y to X.

Regression between microstate metrics/EEG band amplitudes and motor output performance
The second aim of the data analysis was to verify whether a combination of microstate metrics and EEG band amplitudes could predict the performance in the execution of the motor task. This was done by a multivariate regression analysis. The independent variables of the regression model (i.e. the input variables) were obtained as follows. For each trial, the movement period from 5 to 8 s was segmented into six non-overlapping subsequent time intervals of 0.5 s duration. Microstate metrics (mean duration, mean occurrence, transition probabilities) were calculated for each of these intervals. For the same intervals, the mean EEG band amplitudes were calculated using the following procedure. Hjorth montage was obtained in a selection of 14 EEG channels covering parietal and sensory-motor (C3, CPz, C4, P3, P4), visual (P9, POz, P10) and frontal (F3, F4, AFz, FCz, F7, F8) cortical areas. Briefly, for each electrode the weighted average of the neighboring electrodes was taken as its reference [41,42]. Each of these signals was back-forward filtered in the theta (4-7.5 Hz), low alpha (8-10 Hz), high alpha (10-13 Hz) and beta (13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25) Hz) bands using a fourth-order Butterworth filter. The amplitude of the bandpass filtered EEG signals was then obtained by applying the Hilbert transform, and the mean amplitude for each band was calculated in each of the six time intervals of 0.5 s. By doing so, for each time interval we obtained: (a) 56 features for the EEG amplitude (i.e. the mean amplitude values for the 4 frequency bands × the 14 channels considered); (b) k × 3 (i.e. k microstate templates × 3 microstate metrics for each template); These features were calculated for all intervals, all trials and all subjects and were used as independent variables, i.e. as the input of the regression model.
The dependent variables (i.e. the output of the regression model) were the motor performance, estimated by the FV i values for the same time intervals for which the microstate metrics and the band amplitude values had been calculated.
Given their nature, the input variables can be highly correlated with each other. When predicting an output from a set of features, such collinearity could lead to overfitting of the data and poor generalizability of the results [43]. Such limitations can be addressed by both linear and non-linear regression frameworks [44,45]. We opted for a PLS regression model, i.e. a linear regression procedure that assigns weights to all independent variables and exploits a dimensionality reduction of the regression problem [29,46]. The PLS regression is based on the assumption that the measured data are generated by a reduced number of hidden variables that are not directly observable (i.e. the measured data are generated by latent variables that are somehow related to the input variables). The dimensionality reduction in the PLS regression model is implemented by finding linear combinations of all the input variables (i.e. by identifying the so-called components of the model) that optimize the performance of the PLS model. The optimal number of model components is defined by the value of a hyperparameter h that can be either imposed to the PLS model or identified by a datadriven procedure. We opted for the second approach, as explained below in the description of the inner loop of the nested cross-validation (nCV) procedure.
A cross-validation (CV) framework was used to estimate the generalizability of the PLS model. CV works by separating the data in folds and estimating the model on all folds except one-fold (leave-one-out CV [47]). The excluded fold is used to test the performance of the model, which is estimated by the correlation between the PLS model outputs and the measured output values. This procedure is iteratively repeated for different folds.
CV was also used to estimate the optimal value of the hyperparameter that must be set a priori in the PLS model. To this aim, two nested iterations were implemented, consisting of an outer loop and an inner loop (nCV, [48]). The outer loop estimates the generalizability of the model at a group level, whereas the inner loop estimates the optimal value of the hyperparameter h, while also permitting to assess the generalizability of the model at individual level. Figure 2 illustrates how we implemented this procedure with our data.

Outer loop
The number of folds corresponds to the number of subjects (n = 18). For the kth fold (i.e. for the kth subject's data, including microstate metrics, EEG band amplitudes and measured FV i values), data were split in training data (all folds except the kth one) and validation data (the excluded kth fold). To estimate the generalizability of the PLS model, it must be run for each training fold. Therefore, the hyperparameter must be set for each fold through the inner loop of the nCV procedure (see below). Once the outcome of all PLS models for the training data was achieved, the outer loop of the nCV procedure verifies the generalizability of the PLS model on the excluded kth fold (testing data). The outer loop of the nCV procedure is iterated for each fold. After this CV, the PLS model can be generalized by setting an average hyperparameter as the average of the hyperparameters obtained for all subject folds in the inner loops.

Inner loop
To identify the hyperparameter h needed to run the PLS model for each kth fold of the outer loop, only the corresponding training folds are used (i.e. all folds except the kth fold). The training folds are split into two 'sub-folds' , one containing all training folds except the jth fold (these folds compose the training sub-fold) and another one that contains the excluded jth fold (which is the validation sub-fold and consists of the jth subject data). The PLS model is then iteratively run on the training sub-fold for h varying from 1 to the number of the input variables and tested on the validation sub-fold. Once the PLS model has been tested for all h values, the h value corresponding to the best performance is selected for the jth fold (i.e. the validation sub-fold). The inner loop of the nCV procedure is then iterated for each training fold of the inner loop (i.e. for j = 1, …, 17). At the end of all iterations needed to identify the optimal h for testing the kth fold, 17 values for the hyperparameter h are extracted from the training sub-folds. Among these 17 hyperparameter values we choose the one corresponding to the PLS model with the best performance as the optimal hyperparameter to be used in the outer loop to calculate the predicted output and test the PLS model on the kth fold.

Statistical analysis
All microstate metrics are reported as mean values across subjects ± standard error of the mean. To detect the differences between the microstate metrics (duration, occurrence, global explained variance, directional predominance) at rest and during movement, the Wilcoxon rank sum test was separately applied to each microstate. The significance level was set at 0.05 and a Bonferroni correction for multiple comparisons was applied. To detect a significant prevalence of transitions from microstate X to microstate Y (or from Y to X), a single sample Wilcoxon test of the X↔Y directional predominance was performed, checking the difference of the median of these values from zero: a significant positive (or negative) X↔Y value indicated a prevalence of transitions from X to Y (or from Y to X).
A Pearson's correlation analysis between the measured and the predicted FV i values was used to assess the performance of the cross-validated PLS regression model. The Pearson's correlation coefficient r and the null hypothesis significance (p value) were reported. The Pearson's correlation between the measured and the predicted FV i values was repeated for each subject to quantify the performance of the PLS model at single subject level.
The PLS model, similarly to other linear regressors, allows a direct a posteriori evaluation of the relevance of each input feature for the output of the model. Given that the PLS model components are linear combinations of all the input features, it is possible to estimate the general weights of each original feature by reversing the dimensionality reduction transformation. Given that the magnitude of these weights is related to the relevance that each original feature has in the regression model, we analyzed the weights associated with each feature (i.e. each independent input variable) to estimate the significance of the microstate metrics and the EEG band amplitudes in predicting the FV i values. Given that the PLS models were trained to predict the FV i values, positive weights indicate that high values of the independent variables are associated with high FV i values, and thus with low motor performance. Conversely, negative weights indicate that high values of the independent variables are associated with low FV i values, and thus a good motor performance. The statistical significance of the weights was assessed through a randomization procedure, estimating a confidence interval of the null hypothesis of each weight being zero (by means of 10 000 iterations of the same algorithm on randomly shuffled FV i values). The significance thresholds for the weights were estimated considering the 2.5 and the 97.5 percentiles of the randomly shuffled PLS weights distributions.

Microstate templates
According to the KL criterion, we found an optimal number of microstate global templates equal to 4 for both rest and movement conditions. TANOVA between the map topographies at rest and during movement did not show any significant differences for microstates B, C and D (p > 0.500 consistently), whereas the templates of microstate A at rest and during movement were significantly different (p = 0.009). The microstate templates at rest showed topographies similar to those obtained in previous studies [22,49]: a left posterior to right frontal orientation for template A; a right posterior to left frontal orientation for template B; a symmetrical distribution between the two hemispheres, with an occipital to prefrontal orientation for template C; a frontocentral activity for template D ( figure 3). The microstate template A during movement presented a posterior to frontal orientation.
With the optimal number of templates, the global explained variance (mean and standard error) was 83.3 ± 2.2% during movement and 84.3 ± 2.1% at rest.

Rest vs movement comparison: microstate metrics and directional predominance
Given that two different global microstate templates A were found for the rest and movement conditions, for this microstate we only report the values of the metrics in the two conditions in figure 4, whereas we compared the microstate metrics at rest and during movement for microstate B, C and D. The Wilcoxon test applied to microstate duration revealed that microstates C and D were significantly longer during movement than at rest (p < 0.001, consistently), whereas no differences were found for the duration of microstates B (p > 0.5, figure 4(a)). Similarly, occurrences of microstates C and D were higher during movement than at rest (p < 0.001), whereas no difference were found for the occurrences of microstates B (p = 0.100, figure 4(b)).
No significant differences were found for microstate global explained variance.
The values of directional predominance obtained for the rest and movement conditions were compared by the Wilcoxon test, which revealed that the transitions from microstate A to microstate D and from microstate B to microstate D occurred more frequently during movement than at rest (p < 0.001). Please note that the fact that microstate A is different for the rest and movement conditions might influence the statistical results.
The single sample Wilcoxon test against the zerovalue showed that only the transitions from A to D, B to C and C to D during movement had a positive value that was significantly different from zero (p < 0.001, consistently, figure 4(c)), and that only the transitions from A to C during movement had a negative value significantly different from zero (p < 0.001).

PLS regression
Given that four global microstate templates were found with the KL criterion for movement condition, for each time interval we obtained a total of 74 independent input variables for the PLS regression model (12 microstate metrics, 6 directional predominance values and 56 EEG band average spectral amplitudes). When checking for correlations among the input variables of the PLS model, we obtained high correlation values in the upper-(lower-) diagonal elements of the correlation matrix across all subjects, in particular across the EEG channels within the same frequency band and across the microstate metrics ( figure 5). Interestingly, no correlations were found between the EEG band average amplitudes and the microstate metrics.
When performing the PLS regression analysis for predicting the FV i values, the number of the input variables was reduced to a smaller set of uncorrelated components. From the 74 original variables, the nCV procedure identified an average of 6 uncorrelated components, i.e. linear combinations of all the input variables, for the general PLS model. The Pearson's correlation was performed between the FV i values estimated by the PLS model and the measured FV i values. We found a significant positive correlation (r = 0.29, p < 0.000 01, figure 6(a)), confirming that a combination of microstate metrics and EEG band amplitudes is indicative of how constant is a force exerted during an isometric contraction supported by visual feedback. The mean Pearson's r values between the estimated and the measured FV i values calculated for each subject was 0.31 ± 0.12 (mean ± standard deviation; minimum: 0.16, maximum: 0.60, p < 0.001 consistently).  Looking at the estimated weights of the independent variables, we found significant negative values for the microstate duration and occurrence of templates A, B and C, for the global explained variance of templates B, and for the directional predominance from A to D, from B to C and from B to D. A significant positive value was found only for the explained variance of template C. We also found statistically significant (p < 0.05) positive values for the EEG band amplitudes at C4 for the theta band, at F4 for the high alpha band and at FCz and CPz for the beta band. Significant negative weights were found at C3 for the theta, low alpha and high alpha bands, at AFz for the low alpha band, at FCz and P4 for the high alpha band, and at POz, P4 and F7 for the beta band.

Discussion
To elucidate the brain dynamics during the maintenance of a visually driven constant force, we have proposed a multivariate regression model where the collinearity problem was addressed by a data-driven approach. The results obtained confirmed that the proposed PLS model, cross-validated with a nested approach [44,50], reached a good generalizability and succeeded to show the predictive value of EEG microstates and spectral features in inferring the stability of the force exerted during a hand isometric contraction supported by visual feedback. As we discuss in more detail below, our results also suggest that longer duration and higher occurrence of microstates that could be associated with visual and executive control networks corresponded to better contraction performances, suggesting that the visual system and the executive control network play a crucial role in visuo-motor integration.
The experimental protocol was conceived to allow us to select and analyze the time intervals during which only the phenomena related to the maintenance of a stable contraction were present, thus avoiding the interference of confounding factors, such as the response to the 'go' signal, the transient brain activity from rest to movement, and the contraction modulations to reach the requested force level. The functional brain activity was recorded by means of EEG, that provides the high temporal resolution needed to calculate the independent variables that serve as input for the PLS model. The brain dynamics during the performance of the motor task, i.e. the maintenance of a constant force level during an isometric contraction, was represented with a set of different features referring to both local and global brain activity. The local brain dynamics from 14 different sites covering the whole scalp was represented by means of spectral features, calculated as the average EEG rhythm amplitudes at different frequency bands, which are functionally related to the sensorymotor processes. On the other hand, the global brain dynamics was described by means of the metrics of sequences of global microstates. However, both feature types were affected by collinearity, and positive and negative correlations were found among all features. In fact, due to volume conduction phenomena, the electrical field generated by the neuronal assemblies in the cortex is spatially widely spread at the sensor level [51], giving rise to possible correlations among the EEG sensors. This phenomenon, together with possible inter-individual correlations between the EEG amplitudes at different frequency bands, produces a high collinearity among the features, as shown by the elements of the correlation matrix calculated across subjects (see figure 5).
The collinearity among the independent variables used as input features to the PLS model represents a major limitation in building prediction models based on multivariate methods. Furthermore, collinearity can lead to unreliable results because the weights of the model variables can be far from their optimal values and the model may suffer from poor generalizability [45]. To overcome this problem, we applied the PLS regression model within an nCV framework. Starting from the 74 independent variables, we found that 6 independent linear combinations of the microstate metrics and spectral features were optimal for inferring the motor performance during isometric contractions, with a good generalization of the PLS model at both group and individual levels. After having identified the best value of the PLS hyperparameter h in the inner loop of the nCV, it was possible to estimate the general weights associated with each input feature in the generalized PLS model. The value of each weight was then associated with the relevance of the corresponding feature in the PLS regression model, with its sign providing information on the regression direction.
Looking at the sign of the PLS weights, our results showed a negative correlation between the duration and occurrence of microstates A, B and C and the measured FV i values. This means that longer duration and higher occurrence of these microstates were related to low FV i values, i.e. to a better contraction performance. The duration metric has been previously interpreted as being associated with the stability of the neural assembly activities underlying a microstate [49]. Therefore, a longer duration of given microstates associated with a better performance indicates a higher stability of the brain dynamics during the maintenance of a stable force. Similarly, the occurrence of a microstate has been associated with the tendency of the underlying neural generators to become activated [22]. A higher occurrence of a microstate during a good contraction performance could indicate that the general brain activations represented by that microstate underpin the achievement of a good performance. The association of a long duration and a high occurrence with a good motor performance is more evident for microstates A, B and C.
Previous studies attempted to associate single EEG microstate topographies with sets of active brain areas at rest, evidenced by functional MRI [52]. However, no univocal results were achieved [52][53][54]. Other studies localized the sources of broadband EEG scalp topographies [23]. In the light of the results of these studies, microstate B has been associated with the visual area activity [53] and with the visual processing [55]. Therefore, our results on microstate B could be associated with an involvement of the visual system for the visuo-motor integration during task execution.
Moreover, microstate C has been associated with the activity of brain regions belonging to the cinguloopercular network (CON), a well-studied functional system involved in attention control (or executive control), i.e. in the maintenance of behaviorallyrelevant task parameters, of alertness during task performance [56][57][58], and of cognitive resources available for task performance [59][60][61][62]. Given its executive control function, the CON is necessary for volitional motor actions and dexterity, and has been associated with the generation of force and with the motor performance in a visuo-motor tracking task, so that its damage in stroke patients significantly contributes to their motor deficits [63]. Based on these results, higher values of microstate C metrics could be associated with an involvement of brain structures devoted to the executive control function to achieve a good performance in maintaining a constant force level.
The microstate A during movement is a taskspecific microstate. Indeed, its topography, different from the topography of microstate A at rest and with a posterior to frontal orientation, was similar to the topography of the microstate E obtained by Pirondini et al [64] during the holding phase of a reachingand-grasping task. Indeed, these authors found that, in a condition similar to our isometric contraction, the microstate repertoire was characterized by four states: three similar to the microstates found in the resting-state condition (B, C and D) and one microstate, named E, that was task specific. This microstate, as our microstate A during movement, seems to be task-specific and to be related to general motor control strategies during movement execution.
Finally, in previous studies microstate D has been associated with the activity of the right superior and middle frontal gyri and with the right superior and inferior parietal lobules [23], and has been related to the activation of the dorsal attention network [53]. Although the duration and occurrence of microstate D did not show any relevance in the regression model, nonetheless we observed that a higher probability to transit from microstate A and B (associated with the sensory networks) to microstates C and D (associated with higher cognitive control networks) positively correlated with a better performance, as evidenced by the negative weights of A↔D, B↔D and B↔C directional predominance. This result could support an involvement of microstate D in the achievement of a good motor performance.
A direct localization of microstate dynamics would be very interesting and could confirm or weaken the associations made between given microstates and specific brain network activations. However, the localization of brain activity generating EEG microstates is not a trivial problem: the volume conduction effects and the fact that different brain current density distributions can result in the same scalp potential topography makes the solution ambiguous. We used the results of previous studies [23,65] to infer a possible link between microstate dynamics and specific brain networks. It will be interesting to explore the possibility of a direct localization of the observed microstate dynamics in future studies.
The comparison of the microstate metrics between the rest and movement conditions showed a larger representation of microstates C and D during movement (i.e. longer duration and higher occurrence of microstates C and D, and a higher probability to transit from microstate A to D, and from microstate B to C). Although in our case microstate A has a different topography at rest and during movementwhich could influence the statistical comparisonour result is in agreement with a previous finding of Pirondini et al [64] who, by investigating the EEG microstates during the motor control of reachingand-grasping movements, found that the occurrence of microstates C and D significantly increased during the holding phase of the grasping movement. They also showed that, during the execution of voluntary movements, the brain activity evolves passing through a set of spatiotemporal states composed of topographical patterns with spatiotemporal structures similar to those found at rest. Our results also showed similar microstate topographies at rest and during movement (microstates B, C and D), thus confirming that spontaneous brain activity at rest is shared with the basic neural mechanisms underlying motor control.
When we analyzed the sign of the weights of the spectral features in the general PLS model, our results showed a positive correlation between the theta band amplitude over the sensorimotor cortex contralateral to the movement (C3) and a good performance. Accordingly, previous studies outlined a possible functional role of the theta activity in the motor control of the hand. The activation of the theta band has been related to motor plan update, motor learning and neural representations of hand kinematics [66][67][68], situations in which a sensorymotor integration for movement control is required. Moreover, long-range theta coupling between the primary motor cortex and multiple brain regions has been described, also mediated by mechanisms of phase-amplitude coupling with higher frequencies [67,69].
We also found that both the low and high alpha band activities in the contralateral sensory-motor areas were positively correlated with a good motor performance. An ERD of the alpha rhythm was observed during volitional movement [16] and previous studies reported that the alpha band activity is localized in the post-central sensory areas [70,71], and that its modulation during movement is related to the processing of sensory information [72,73]. In a visuo-motor task similar to the one used in our study, Porcaro et al [74] observed that a significant higher alpha ERD occurred in trials with high performance but only during the first second of the movement, when a huge sensory processing is needed to adjust the level of the exerted force to the required force level. In our case, during the intervals with a bad motor performance, a huge visual mismatch between the level of the exerted force and the required force level was present, resulting in a higher demand for sensory processing to adjust the motor output. This higher demand was related to a lower alpha amplitude.
Interestingly, while we observed that the EEG band amplitudes at C3 (i.e. the brain activity over the motor cortex contralateral to the movement) were associated with a good performance (i.e. with a stable exerted force), the EEG band amplitudes at C4 and F4 (i.e. over the ipsi-lateral motor/premotor areas) were associated with a higher variability of the exerted force, hence with a poor performance. This result might be compatible with the 'interhemispheric competition' model [75]. According to this model, the contralateral hemisphere controls the motor execution while suppressing the activity of the ipsilateral hemisphere to reduce putative interference of ipsilateral descending pathways believed to degrade motor performance. Within this framework, our results during movement seem to indicate the onset of a kind of inter-hemispheric competition aiming at achieving an optimal motor performance.
Differently from other studies, where a positive relationship between the beta band activity over the primary motor cortex and a stable motor output during a visuo-motor task was found [12], we observed that the values of the weights in the beta band at C3 were not significantly different from zero. However, the beta band amplitude over FCz, a site corresponding to the fronto-central brain areas, was negatively correlated with the motor performance, confirming the role of the beta activity in maintaining the current sensorimotor state [10]. It is known that fronto-central brain areas, as the supplementary motor area and the anterior cingulated cortex, are generally involved in self-generated action [76].
We also found positive and negative correlations of the alpha and beta band amplitudes with motor performance in parietal sites, thus suggesting that the achievement of a good motor performance is also supported by the involvement of different brain areas functionally connected for visuo-motor integration. In fact, previous transcranial magnetic stimulation studies [77] documented functional connections between contralateral primary motor cortex and the parietal cortex and its sub-regions (such as posterior parietal cortex and different portions of the intraparietal sulcus). As well, transcranial direct current stimulation of the right parietal cortex (P4) was shown to induce changes of the primary motor cortex excitability during motor imagery and action observation [78].
Interestingly, we did not find a direct correlation between the microstate metrics and the spectral amplitudes, suggesting that both types of features contribute, although differently, to describe the brain dynamics that support the achievement of a good motor performance during the maintenance of a stable visually guided force. Given that microstate metrics and EEG amplitudes characterize brain dynamics at different scales, the absence of correlation among these features supports the concept that they represent complementary aspects of brain processing. It is worth noting that, based on the results of the nCV of the PLS model, microstate metrics and brain rhythm amplitudes together seem to explain the brain dynamics underpinning the maintenance of a visually guided stable motor output not only at a group level, but also at an individual level. This finding could play a crucial role for a better understanding of the motor system in single trials or in realtime applications, having an impact on the progress of brain computer interface systems [79], in the study of motor control development [80] or of motor control changes due to aging [81].

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.