Blindly separated spontaneous network-level oscillations predict corticospinal excitability

Objective. The corticospinal responses of the motor network to transcranial magnetic stimulation (TMS) are highly variable. While often regarded as noise, this variability provides a way of probing dynamic brain states related to excitability. We aimed to uncover spontaneously occurring cortical states that alter corticospinal excitability. Approach. Electroencephalography (EEG) recorded during TMS registers fast neural dynamics—unfortunately, at the cost of anatomical precision. We employed analytic Common Spatial Patterns technique to derive excitability-related cortical activity from pre-TMS EEG signals while overcoming spatial specificity issues. Main results. High corticospinal excitability was predicted by alpha-band activity, localized adjacent to the stimulated left motor cortex, and suggesting a travelling wave-like phenomenon towards frontal regions. Low excitability was predicted by alpha-band activity localized in the medial parietal–occipital and frontal cortical regions. Significance. We established a data-driven approach for uncovering network-level neural activity that modulates TMS effects. It requires no prior anatomical assumptions, while being physiologically interpretable, and can be employed in both exploratory investigation and brain state-dependent stimulation.


Introduction
Transcranial magnetic stimulation (TMS) applied to the human neocortex produces highly variable effects [1][2][3].This variability can be partially attributed to the dynamic nature of neural activity in the stimulated brain area.Rather than treating the TMS readout as a true effect obscured by noise, this can be seen as the result of the interplay between the stimulation effect and the brain's endogenous neuronal activity [4].Such a conceptual approach to TMS effects allows exploiting variability in the readout to study dynamic brain states [5].In this study, we explored this approach in combination with machine learning (ML) analysis techniques to identify functionally relevant patterns of cortical activity in relation to TMS effects.
Previous studies within the sensorimotor network have associated fluctuations in TMS effects with the state of oscillatory cortical activity recorded with electroencephalography (EEG).The power and phase characteristics of neuronal oscillations recorded just before or at the onset of stimulation have been linked to the level of TMS-induced excitation of corticospinal pathways, represented by motor evoked potentials (MEP) [6].The effects of endogenous neuronal oscillations on TMS-evoked activity may originate in the dynamics of local and global cortical excitability, i.e. the probability that a given neuronal population will respond to an input signal [7].Furthermore, excitatory and inhibitory connections from other regions of the functional network in study may affect the excitability of the stimulated area.The connectivity state between communicating cortical areas may be inferred from their oscillatory signals by measuring the alignment of their phases [8].Indeed, in primates, the spiking rate of individual neurons was found to depend on the phase-coupling between local field potentials in distal regions [9].Thus, relative phases across brain regions, in addition to local oscillatory power and phase, may partially explain variability in the MEP amplitude.The variability of TMS effects offers a unique opportunity to investigate state dynamics of the brain network activity.
Despite active ongoing research, neuronal processes modulating TMS effects within the sensorimotor network are still largely unknown.In the context of EEG-TMS studies, some of the obstacles on this path are low signal-to-noise ratio of EEG oscillations, mixing of source signals in scalp recordings due to volume conduction, and scarcity of prior knowledge about functionally relevant cortical sources.To overcome the first two issues, source activity can be reconstructed through spatial filtering of sensor signals, wherein a weighted average of sensor signals is taken [10].However, commonly used 'model-based' source reconstruction techniques require prior assumptions about source locations as well as a forward model (i.e. the signal's sourceto-sensor mixing process) [11].In the absence of either, spatial filtering can be achieved with 'datadriven' Blind Source Separation techniques, such as Common Spatial Patterns (CSP).CSP is designed to separate a multivariate EEG signal into components that are most distinguishing between discrete experimental conditions or outcomes [12,13].CSP is commonly used in the field of brain-computer interfaces to decode right-vs.left-hand motor imagery from periods of spatially specific event-related desynchronization detected in EEG signals [14].In our study, we used a variant of CSP, called analytic CSP (aCSP), that is particularly suitable for studying oscillatory phenomena as the analytic signal (derived using the Hilbert transform) simultaneously encodes instantaneous phase and amplitude using complex numbers, and this allows a separation of multivariate signals based on both amplitude relationships and also phase relationships.Signal components produced by aCSP, commonly referred to as spatial patterns, can capture dynamics of both local oscillatory amplitude (e.g.standing waves) as well as network-level phasespecific communication between neuronal populations (e.g.phase-coupling phenomena and travelling waves).Employing these components as features in an ML classifier quantifies their relevance for the experimental conditions.In principle, aCSP provides an opportunity to uncover neuronal correlates of targeted brain function, while requiring little prior knowledge and/or assumptions about the nature of the neuronal activity in question.
Within the domain of EEG-TMS, aCSP can be employed to detect neuronal processes, which interact with the effects of TMS and thus predict stimulation outcomes.Previous studies have shown applicability of ML methods to EEG-TMS-based brain state identification [15,16].We applied aCSP to an EEG-TMS dataset obtained in the course of single-pulse TMS of the primary motor (M1) cortex.We used EEG components extracted with aCSP as features in an ML classifier to predict MEP amplitudes from pre-stimulus EEG signals.Furthermore, we examined the components to identify functionally relevant cortical activity and describe its spatial, temporal, and spectral characteristics [17].Overall, this paper proposes a new data-driven approach to studying the variability of TMS effects and their relationship with brain activity that could lead to a better understanding of the underlying neuronal processes.Such an approach is applicable at an exploratory stage of an investigation, as well as within the domain of brain state-dependent stimulation both in research and clinical application.

Dataset
The dataset used in this study consisted of 20 EEG-TMS recordings.Experiments were performed on right-handed healthy adult participants with no known neurological conditions (12 females, 8 males, mean age ±SD = 26 ± 4).All participants provided written informed consent prior to participation.The study was approved by the ethics committee at the Faculty of Medicine in the University of Tübingen (approval ID: 810/2021BO2) and conducted in accordance with the Declaration of Helsinki.The EEG-TMS recordings were acquired previously for other purposes with slight variations in the experimental protocol.Single TMS pulses (1000-1200 pulses) were applied over the hand knob area in the left M1 at 110% of the resting motor threshold (RMT) at 2-3 s intervals with random jitter (2 ± 0.25, 2.1 ± 0.1, 3 ± 0.5 s, depending on the protocol of the given recording).EEG with 126 channels (positioned following the international 10/10 placement system) recorded continuous signal from the scalp with a 5 kHz sampling rate, while two bipolar EMG channels recorded activity from the abductor pollicis brevis (APB) and the first dorsal interosseous (FDI) muscles of the right hand.The experimental procedure, the acquired dataset, as well as the data preprocessing pipeline are described in more detail elsewhere [15] (the data acquisition description for participants 1-9 is described in Zrenner et al [18]).

EEG preprocessing
EEG data was preprocessed using the methods described in Metsomaa et al [15].Briefly, trials were epoched around the TMS pulse, raw EEG signal within a 1.5 s window (1.500-0.005s) before each TMS pulse was downsampled to 1 kHz, after which slow trends were removed from the signal.Noisy channels and trials were excluded from the data based on their deviation from the respective median noise level.Then, eye movement artefacts were removed from the signal with Independent Component Analysis.Cleaned EEG data were filtered in the 8-30 Hz band with a 6th-order Butterworth filter and downsampled to 250 Hz.The frequency range of the bandpass filter spanned the alpha-and betarange, both of which are known to be functionally relevant oscillatory frequency bands within the sensorimotor network [19].Finally, the signals were transformed into their analytic representation with the Hilbert transform.The signal within a 0.5 s window (0.505-0.005 s) preceding the TMS pulse was used for the subsequent analysis, leaving 124 time samples in each epoch.The time window was selected such that it included at least four cycles of each retained frequency [20].

EMG preprocessing
EMG signals from the APB and FDI muscles of the hand were preprocessed in the following way: continuous EMG signals were separated into 1 s-long epochs centered at each TMS pulse.Epochs containing pre-innervation within the 300-ms pre-stimulus window of the EMG (defined as a maximum peakto-peak signal exceeding an individual threshold set between 30 and 40 µV) were excluded from further analysis.Slow trends and TMS-related artefacts were removed from the remaining epochs (for a detailed description of the preprocessing procedure, see Metsomaa et al [15]).MEP amplitudes were calculated on the clean EMG signals as a peak-to-peak amplitude distance in the 18-55 ms window after the TMS pulse.MEP amplitudes from the muscle with the higher average amplitude value for a given subject were selected for further analysis.All trials were divided into a 'High' and a 'Low' corticospinal excitability condition (henceforth referred to simply as 'High' and 'Low') based on the respective MEP value.In order to do so, while taking into account possible slow trends in the amplitudes across the experimental session, a dynamic baseline was defined as a moving median of 150th order across successive trials.The MEP values were labelled as 'High' or 'Low' depending on whether they were above or below their respective baseline.

aCSP
aCSP was applied to EEG data, following the approach described in Falzon et al [13] (figure 1).aCSP decomposes multivariate data into a set of components using generalized eigenvalue decomposition (GED).This method takes complex spatial covariance matrices for each condition as input and generates a set of eigenvectors and eigenvalues.The eigenvectors are used as spatial filters to extract signal components that account for the maximum variance in one condition and the minimum variance in the other.These components can be considered as reconstructed source-level neuronal activity that exhibits the difference between the experimental conditions.aCSP analysis was applied to each subject's EEG data in the following way.
All available trials in a given condition were ranked based on their respective MEP values (see EMG preprocessing).200 trials with the highest MEP values in the 'High' condition and the same number of trials with the lowest MEP values in the 'Low' condition were selected, leaving 400 trials for the subsequent analysis.This selection aimed to maximize the separability of the two conditions.
Within each condition, spatial complex-valued covariance matrices were calculated from the prestimulus analytic EEG signals and averaged across trials.From each trial's EEG epoch contained in an n × m complex-valued matrix denoted as X, with n being the number of EEG channels and m being the number of time samples, a normalized n × n covariance matrix R was calculated as: R = XX * tr(XX * ) .The denominator in the equation is the trace of the covariance matrix or the sum of the squares of the samples from each channel and X * denotes complex conjugate transpose of X.The two averaged covariance matrices were then used for the aCSP analysis.
aCSP was performed by means of GED for each condition separately, maximizing signal variance in the chosen condition while minimizing the total variance (represented by a sum of the two averaged covariance matrices from both conditions).To prevent overfitting to noise, the GED was regularized with a coefficient weighted by channel-wise variances.The regularization coefficient was selected via a cross-validation (CV) procedure as a value between 1e −8 and 1e −1 (see Classification and CV), and the channel-wise variances were averaged across all analyzed trials from both conditions pulled together.
The output of aCSP is represented by an n × n matrix, where each column is an eigenvector, accompanied by a set of n corresponding eigenvalues.These eigenvectors serve as spatial filters for the sensor-level EEG signals.The eigenvalues represent the proportional differences in the amount of variance explained by each eigenvector between the two conditions.We selected eigenvectors with the largest eigenvalues, maximizing signal variance for one of the two conditions.Between 2 and 6 spatial filters were chosen for further analysis, with an equal number of filters for each condition.The number of filters was selected via the CV procedure and differed across the CV folds (see Classification and CV).Since input covariance matrices were complex-valued, the resulting aCSP filters were also complex-valued.
The variance of each aCSP component was then used as a predictor for the excitability condition label.Specifically, the variance of each component in a given trial quantifies the power of that component within the EEG signal in that trial.The given component's predictor feature p was calculated from the single spatial filter contained in a n × 1 complex-valued vector f as: p = |f * Rf |.This yielded one feature value per trial for each aCSP component (i.e.2-6 values per trial, depending on the number of filters in a given CV fold).
We performed a few variations of the analysis in order to characterize the predictive components.In order to verify, whether the predictive component is time-locked to the stimulation event, we repeated the analysis with both variance and phase of the component time courses as predictor features in LDA.While variance of a time course served to quantify presence of the component in the analyzed time window, phase of a time course served to quantify the extent to which the activity of the predictive component was time-locked to the stimulation event, i.e. whether the time courses were aligned across trials in their phases with respect to the stimulation event.The phase was derived from the spatially filtered complexvalued time course at the last time sample of the analyzed window (12 ms before the pulse) and was transformed into sine and cosine of the angle of the complex value before being passed on to LDA.Next, we performed the analysis with only phase features as predictors in LDA.In both cases sine and cosine of the angle were passed to LDA as two separate features.Finally, we tested whether phase-shifted network activity played a role in the prediction of the excitability state.Instead of using analytic signals, we performed CSP on real-valued EEG signals that did not undergo Hilbert transform.Real CSP extracts purely instantaneous activity (i.e.changes in signal variance happening instantaneously across the scalp), while aCSP extracts both instantaneous and phaseshifted activity (i.e.changes in signal variance phaseshifted across the scalp).In addition, we compared the power of real and aCSP components within the same trials for each subject by calculating the filters on the same training set (320 trials), applying them to the same test set (80 trials), calculating logarithm of the variance of each component within each trial and performing Pearson correlation between CSP and aCSP components across test trials.

Classification and cross-validation
A classification analysis was performed to test whether the obtained aCSP components were predictive of the excitability condition.In each iteration of the analysis, 400 trials from a single subject's data were randomly divided into a training and a test set in a 4:1 ratio, i.e. with 320 trials in a training set and 80 trials in a test set.Both the training and the test sets included an equal number of trials from each condition.The aCSP filters were generated using the training set, and they were then applied to both the training and the test sets to create predictors for the classification.Regularized Linear Discriminant Analysis (LDA) with automatic hyperparameter optimization was then trained on the training set and applied to the test set to predict the condition labels from the variance of the aCSP components.The percentage of correctly classified test set trials was taken as the classification accuracy.
To ensure that the classification results were not driven by possible outliers in the randomly assigned test set, the results were calculated on and averaged across various training-test partitions of trials.This was implemented via a 5-time 5-fold CV procedure (figure 1).Each subject's data were randomly split into five equal folds, with one fold assigned as a test set.The whole analysis including aCSP and LDA was repeated with different training-test set re-assignments until all available folds had been used as a test set once (i.e.five times).Furthermore, the partitioning of data into folds was repeated five times to average out the effects of randomness in the data splitting process.The average classification accuracy across all 5 × 5 repetitions of analysis was taken as the classification accuracy of the given subject.
Furthermore, the CV operated on two levels.The outer layer was dedicated to the estimation of a given subject's classification accuracy and was performed as described above.The inner layer of CV was dedicated to the selection of analysis hyperparameters, i.e. the number of aCSP filters (2, 4 or 6 filters) and the value of the aCSP regularization coefficient (1e −8 , 1e −6 , 1e −4 , 1e −2 , and 1e −1 ).The hyperparameters were selected via a 5-fold CV procedure, which was performed anew for each fold iteration of the outer layer.
The combination of hyperparameters that yielded the highest average classification accuracy across the five folds of the inner layer was used on the outer layer.In this way, the optimal hyperparameters were estimated individually for each CV fold.

Spatial patterns analysis
The aCSP filters can be viewed as inverse operators to retrieve neuronal source activity from multidimensional EEG signals.The filters can be transformed into spatial EEG patterns, also known as activation patterns or topographies, which are then equivalent to forward models.These patterns reflect how the source signal projects onto the sensor space (i.e.source-tosensor spatial mapping), and are in principle neurophysiologically interpretable.
The spatial patterns could not be obtained directly from the classification analysis due to the use of a CV procedure with varying hyperparameters and subsets of data.Consequently, the spatial patterns used for interpretation were obtained in a separate analysis procedure and, therefore, were associated with but did not directly correspond to either the aCSP components used in the classification analysis or the classification results.The spatial patterns were obtained in the following way.To obtain the spatial patterns for a given subject, the aCSP filters were calculated on all 400 trials, without separation into training and test sets.The regularization coefficient value was obtained from a CV-fold of the main analysis with the highest classification accuracy.The spatial patterns were calculated using the method described in Haufe et al [17].The complex-valued n × n matrix W, with columns being the filters, was transformed into a complex-valued n × n matrix A, with columns being spatial patterns, as: A = W −⊤ , where −⊤ denotes the complex conjugate transpose of the inverse.Before averaging across subjects, each pattern was individually normalised to unit-norm.
We selected a single pattern with the largest eigenvalue, i.e. the most distinct between the two conditions, for each condition and each subject for all further 'interpretation' analysis.The patterns were calculated only for the subjects that yielded classification accuracy above chance level (see Statistical analysis).Since both the filters and the associated spatial patterns were complex-valued, the magnitude and phase parts of the pattern topographies were visualized separately [13].The magnitude part represents the distribution of amplitude across sensors, and was derived by taking the absolute values of the complex-valued pattern.The phase part represents the distribution of phase shifts across sensors.While the complex-valued aCSP filters used in the main analysis, as well as the spatial magnitude patterns derived from them, were reference-free (i.e. had an averaged reference), the phase patterns derived at the interpretation stage were re-referenced to particular channels in order to visualize relative phase shifts across the scalp.Before deriving the phase pattern, the complex-valued 'High' pattern was referenced to the FCC3h channel, so that the value at each channel quantifies the phase shift between that channel and FCC3h.Meanwhile, the 'Low' pattern was referenced to the Cz channel.The phase pattern was then derived by computing the angles of the complexvalued pattern.The choice of a reference channel was to some extent arbitrary, selected such that they result in visually smooth phase shift distributions across the scalp in the respective group-averaged pattern to facilitate visual interpretability.All missing channels were interpolated with spherical spline interpolation before computing magnitude and phase patterns.

Time and frequency analysis
Filtering EEG signals with an aCSP filter produces a time course of the extracted signal component.This time course can be assessed in terms of its temporal dynamics and spectral characteristics.The following analysis was performed only on the subjects with significant classification results.A 1.25 s-long prestimulus analytic EEG signal with a 250 Hz sampling rate was divided into a training and a test set (160 and 40 trials, respectively, within each condition).The aCSP filters were then calculated on a 0.5 s window of the training set and applied to a 1.25 s window of the test set.The power spectra of the time courses were estimated within an 8-30 Hz frequency range, averaged across trials, and then across subjects.To visualize the temporal dynamics of the components, the envelope of the filtered signals, reflecting the amplitude as a function of time, was computed by taking the absolute value of the analytic signal and averaged across trials of the test set, and then across subjects.Before averaging across subjects, the envelope was normalized to unit-norm within each component but across trial-average time courses from both conditions pulled together, for visualization purposes.
To further explore the spectral specificity of the discriminative components, we repeated the main analysis on signals that were bandpass-filtered to different frequency bands: 4-8, 8-13, 13-30, and 30-40 Hz within the same time window of 0.5 s preceding the TMS pulse.The frequency bands were not equal in size but rather were chosen to correspond to theta, alpha, beta, and low gamma-rhythms, respectively.Of note, theta-and low gamma-bands were outside of the frequency spectrum of the main analysis but were still included for comparison.To explore the time specificity of the components, the main analysis was repeated at different latencies of the pre-stimulus signal, overlapping by 250, 1250-750, 1000-500, and 750-250 ms before the TMS pulse.The signals were bandpass-filtered to 8-30 Hz, as in the main analysis.

Statistical analysis
To account for possible deviations of the data from normality, the significance threshold was determined by performing a permutation test on each of the 20 subjects.For that, the aCSP procedure was repeated with the following modification.After the trials were selected based on their MEP amplitudes (see aCSP) but before proceeding with the aCSP analysis, the condition labels were randomly permuted across the selected trials.The permutation procedure and the subsequent analysis were repeated 1000 times, and the resulting classification accuracy values formed a null distribution of the classification results.The accuracy at the 95th percentile of the distribution marked an upper confidence limit for a given subject and its average value across subjects was taken as a significance threshold for all subjects.Of note, the permutation test was performed for the main analysis (i.e. with variance of aCSP components used as a predicting feature) but not for any further analysis that used other predicting features, for which statistical significance was not evaluated.
Since aCSP analysis was performed on each subject independently, there is no imposition on the extracted signal components to represent the same neurophysiological phenomenon across subjects in terms of function or spatial localization.The only explicit commonality between them is the predictiveness over the excitability condition.However, the spatial similarity between aCSP patterns from different subjects would indicate the physiological validity and generalizability of the individually derived patterns on a population level.The spatial similarity across the aCSP patterns was statistically tested with a correlation analysis on a subset of subjects with statistically significant classification results.For each selected pattern, channels excluded at the preprocessing stage were interpolated before the analysis.Then, correlation coefficients were calculated between each individual magnitude topography and the group-average topography.The average correlation coefficient was then taken as a measure of similarity.The analysis was performed separately for each excitability condition.
The statistical significance of the result was evaluated with two permutation tests, by shuffling either EEG channels within the selected patterns or the selection of the patterns as such.The first test compared the similarity among the patterns against the similarity between random spatially uncorrelated sets of values.Before the correlation analysis, channels were randomly permuted within each topographical map.The procedure was repeated 10 000 times, and the correlation value at the 95th percentile of the resultant null distribution was taken as a significance threshold.The second test compared the similarity among the patterns against other aCSP patterns, thereby testing the uniqueness and salience of the patterns associated with the highest eigenvalues.
When GED is performed, the same number of filters is generated for each condition as there are channels present in the EEG data, although only the ones with the highest eigenvalues are utilized for most analysis.Here, instead of permuting channels within a tested pattern, the remaining aCSP patterns from the same condition were used.On each iteration, a pattern was randomly drawn from a complete set of patterns within a given condition from a given subject.Correlation coefficients were calculated between randomly drawn patterns and their average pattern, and finally a mean correlation coefficient was taken.This procedure was repeated 10 000 times, and the correlation coefficient value at the 95th percentile of the obtained distribution was taken as a significance threshold.

Validation analysis
In order to validate the results, we applied the same analysis to a different EEG-TMS dataset.This dataset consisted of recordings from 10 healthy righthanded adult participants who did not have any known neurological conditions.Single-pulse TMS was applied over the hand knob area in the left M1 using 800 pulses at an intensity of 110% of RMT, with an inter-stimulus interval of 2.25 s and a random jitter of ±0.125 s.EEG was recorded throughout the experiment by a 64-channel system with a 5-kHz sampling rate, while EMG activity from the APB and FDI muscles of the right hand was measured with two bipolar channels.The preprocessing and analysis pipeline used in the main analysis was applied to this dataset the same way as described above for the primary dataset.

Example of an individual analysis pipeline
For illustrative purposes, in the following section we will present the analysis steps and results using an individual subject as an example case (figure 2).The preprocessed EEG epochs were categorized into either the 'High' or 'Low' condition based on their corresponding MEP amplitudes (figure 2(A)).From each condition, 200 trials were selected for further analysis.The trials with the highest and lowest MEP amplitudes above or below the moving median, respectively, were chosen.aCSP was then performed on the selected data, generating spatial filters targeted at each condition (figure 2(B)).The EEG data was subsequently spatially filtered to isolate the signal components, and the variance of these components in each trial served as predictors of the condition label in LDA classification.We employed the nested CV procedure, repeating the aCSP + LDA analysis multiple times with different subsets of trials and different hyperparameter values.The overall prediction accuracy for each subject was determined by calculating the average classification accuracy across all CV rounds.
For the interpretation stage, all 400 trials of EEG data were used in the GED, generating signal components without conducting any classification analysis.We selected only the components with the highest eigenvalues (i.e.most discriminating between the excitability conditions) for further interpretation.The spatial, temporal and spectral characteristics of these components are visualized in figures 2(B)-(F).As the filters were complex-valued, the magnitude (i.e.absolute value) and phase topographies are visualized separately (figure 2(B)).Although the filters themselves are not directly visually interpretable, they can be transformed into spatial patterns that allow for physiological interpretation.The patterns are visualized as pairs of magnitude and phase maps (figure 2(C)).When considered together, they describe the progression of an oscillatory signal across the scalp.The dynamic nature of the complexvalued spatial pattern can be alternatively depicted as the change in voltage distribution as a function of phase (figure 2(D)).Multiplying the pattern with a generic unit-amplitude oscillation projected from the source offers a different view of the voltage dynamics.In further text, we will restrict pattern visualization to magnitude and phase topographies, as shown in figure 2(C).
Finally, we spatially filtered EEG signals with the aCSP filters to derive time courses of the associated signal components.The amplitude of these time courses corresponds to the component's presence in the EEG signal at a given point in time (figure 2(E)).Importantly, the phase of the filtered signal at any given latency does not play a role in the MEP prediction, only the variance of the component across the analyzed pre-stimulus window does.The filters were created based on 0.5 s-long signals but were applied to longer 1.25 s epochs from a set of trials not used when generating the filters.For this subject, the presence of the components in the signal fluctuates with time but does not exhibit any change in the temporal proximity to the pulse.However, the presence of the 'High' component is consistently more evident in the trials that resulted in high MEP amplitude as opposed to the trials that ended with low MEP amplitude.This was not the case to the same extent for the 'Low' component.Additionally, the time courses were decomposed into their spectral representation (figure 2(F)).Here the spectral peak in the alpha-frequency band is evident, for both the 'High' and the 'Low' components.

Classification accuracy
To assess the predictive value of the aCSP components for MEP amplitude, we used the variance of the components in each trial as features in a classification test (see Classification and CV).The accuracy of the classification, measured as the proportion of the correctly classified trials with respect to their predefined excitability labels, served as an indicator of prediction success.The average classification accuracy across 20 subjects was 68% ± 8% (mean ± SD), ranging between 57% and 91%.The statistical significance of the individual classification results was evaluated by establishing confidence limits from a null distribution (see Statistical analysis).The null distribution's median was at 50% for all subjects, while the group-average upper confidence limit was 59% ± 0.5%.With the significance threshold set at 59%, the excitability condition was successfully predicted for 19 out of 20 subjects.This accuracy was achieved when the number of used aCSP components was allowed to vary across CV folds between 1, 2, or 3 highest components per each experimental condition (2, 4, or 6 components in total).When restricting the analysis to the single highest aCSP component per condition (2 components in total), the average classification accuracy was 66% ± 9%, with 16 out of 20 subjects exhibiting significant prediction accuracy.
We repeated the analysis with the real-valued EEG signals, rather than their analytic representation, which would correspond to the standard CSP approach.Thus, we implicitly tested whether addition of phase-shifted network activity to instantaneous activity improved prediction of the excitability state.The average classification accuracy across 20 subjects was 68% ± 8%, ranging between 55% and 92%, making the classification success identical between CSP and aCSP.We further repeated the aCSP analysis, employing both the component's variance across the trial as well as the instantaneous phase at the end of the analyzed time window (12 ms before the pulse), as predictor features in LDA.Thus, we tested whether the effect of the aCSP component on the state prediction was time-locked to the stimulation event.Since the phase is not a linear measure, it was included in the LDA in the form of sine and cosine of the phase as two separate features.The average classification accuracy across 20 subjects was 67% ± 8%, ranging between 54% and 92%, making the classification success identical with the main results of aCSP, which included only variance as the predictor.Finally, we ran the aCSP analysis with only the instantaneous phase at the end of the window of the component time course as a predictor in LDA.The average classification accuracy across 20 subjects was 51% ± 2%, ranging between 48% and 58%, making it essentially chance level prediction.

Spatial patterns
To assess the spatial similarity of the aCSP components, we derived spatial magnitude maps from the individual components (see Spatial patterns analysis).We averaged these maps across subjects and measured the spatial correlation between each subject's map and the average map within each condition (see Statistical analysis).To evaluate the statistical significance of correlation, we calculated confidence intervals (CI) using two approaches: random permutation of channels (referred to as CI-channel) and random selection of a pattern from each subject's full set of derived aCSP patterns for a given condition (referred to as CI-pattern).
The magnitude map represents how the component's underlying sources are projected onto the scalp.The average spatial pattern of the 'High' excitability condition exhibited distributed localization in the left central-parietal, left frontal-central, right frontal and occipital areas (figure 3(A)).The analyzed individual patterns exhibited significant spatial correlation, although they were no more correlated with each other than a combination of any other individual patterns generated by the aCSP for the 'High' condition (Pearson's r = 0.3, CI-channel = 0.03, CIpattern = 0.36).Upon visual inspection, a similar pattern recurred in 6 out of 19 subjects, primarily located in the left central-parietal area (figure 3(C)).
In the 'Low' excitability condition, the average spatial pattern was localized in the right parietaloccipital, medial frontal and bilateral temporal areas (figure 3(B)).The analyzed patterns also exhibited statistically significant similarity, but to no greater extent than other patterns from the same condition (Pearson's r = 0.25, CI-channel = 0.03, CIpattern = 0.34).Upon visual inspection, a similar pattern repeated in 6 out of 19 subjects, localized in the medial parietal-occipital and frontal areas (in a different subset of subjects compared to the 'High' condition, figure 3(D)).
To assess the phase patterns of the components, we calculated phase maps for subjects with similar magnitude patterns in 'High' (figure 3(E)) or 'Low' (figure 3(F)) conditions (N = 6 in each subset).The phase values represent the phase shift in each channel with respect to the reference, and their signs indicate the direction of phase progression.The signs are arbitrary and depend on the choice of the reference EEG channel, from which the phases were subtracted.The average phase pattern in the 'High' condition revealed a phase shift relative to the phase in the FCC3h channel along the posterior-anterior direction slanting toward the vertex (figure 3(E)).When re-referenced to more posterior channels (e.g.CP3), the direction of the phase shift changed, suggesting the presence of a travelling wave along the posterioranterior path through the left central area of the stimulated hemisphere.For the 'Low' condition, the topographical phase distribution showed a phase shift along the posterior-anterior direction relative to the phase in the Cz channel, indicating a travelling wave along the posterior-anterior path through the midcentral area (figure 3(F)).
Since real and aCSP components performed equally well in the MEP amplitude prediction, we verified whether CSP and aCSP isolated the same EEG components.We did so by running one iteration of both CSP and aCSP analysis on the same training set (320 trials) and applying the first filters from both of them to the same testing set (80 trials).Then we calculated Pearson correlation coefficient between the log of variance of CSP and aCSP components across test trials in each subject.The components were significantly correlated in all subjects in both conditions (p < 0.05), with an average correlation coefficient across subjects in the High component of 0.79 ± 0.2, while the average in the Low component was 0.81 ± 0.23.

Time and frequency analysis
We evaluated the spectral composition of the aCSP components by analyzing the power spectra of the spatially filtered EEG signals (see Time and frequency analysis).To distinguish the spectral characteristics of the components from the intrinsic spectral properties of the non-filtered EEG signals, we applied the same filters to the EEG data from both experimental conditions.Across subjects, the frequencies in the alphafrequency range dominated the spectrum, regardless of the applied filters or the condition, which the EEG signals belonged to (figures 4(A) and (B)).The ). See also figures S2 and S3.(C) Group-average magnitude topography for the 'High' spatial pattern calculated on a subset of subjects with a similar 'High' topography (N = 6).(D) Group-average magnitude topography for the 'Low' pattern calculated on a subset of subjects with a similar 'Low' topography (N = 6).(E) Group-average phase topography for the 'High' pattern calculated on a subset of subjects with a similar 'High' magnitude topography (N = 6, same subset as in (C)).The colors represent the phase shift (in degrees) with respect to a reference channel (FCC3h channel, indicated with a black dot).All phases were subtracted from the phase in the reference channel.The channels exhibiting the magnitude below the median of the magnitude distribution were masked (in grey).(F) Group-average phase topography for the 'Low' pattern calculated on a subset of subjects with a similar 'Low' magnitude topography (N = 6, same subset as in (D)).The colors represent the phase shift (in degrees) with respect to the reference channel (Cz channel).
difference in power of the component between the EEG data from the two conditions was more prominent with the 'High' component than with the 'Low' one.Furthermore, the difference in power between the components was more pronounced with the EEG data from the 'High' rather than 'Low' condition.
To examine the role of oscillatory activity in different frequency bands in the effect of TMS, we conducted the analysis using narrower bandpass-filtered EEG signals (figure 5).The highest prediction accuracy was achieved when the signals were filtered in the beta-frequency band (66% ± 8%).However, this accuracy result is still not as high as with the broadband 8-30 Hz signal used in the main analysis.Analysis on the theta-, alpha-, and low gammafrequency band-filtered analytic signals resulted in accuracies of 60% ± 7%, 62% ± 7%, and 64% ± 9%, respectively (figure 5(A)).We hypothesized that the lower prediction accuracy with narrower spectral filtering may be due to signal distortion.To test this, we further divided the beta-band into two narrower sub-bands, low beta (13-22 Hz) and high beta (22)(23)(24)(25)(26)(27)(28)(29)(30), and repeated the analysis.The classification accuracy decreased to 64% ± 8% for both subbands, which was still 2% higher than the accuracy obtained with the alpha-band signal (which had a still narrower spectral filter), but equivalent to the low gamma-band results (which had the same filter width).The topographical maps appeared to exhibit more focal patterns when narrower frequency bands were isolated compared to the results of a broadband analysis (figure 5(B)).
For consistency, we averaged the topographical maps from each frequency band across all analyzed subjects, regardless of the individual statistical significance of the classification accuracy.The spatial patterns revealed that within the alpha-frequency band, the 'High' condition was predicted by signals from left central-parietal region, while the 'Low' condition was predicted by signals from medial parietal-occipital and frontal locations.Within the beta-frequency band, the 'High' condition was predicted by signals localized in the left frontalcentral, right frontal and occipital areas.The 'Low' beta-component was localized in bilateral temporal regions.Notably, we observed that both 'High' and 'Low' spatial patterns derived from the broadband 8-30 Hz signal (figure 3(A)) appeared to be superpositions of respective alpha-and beta-specific patterns (figure 5(B)).
We examined the temporal dynamics of the aCSP components within the pre-stimulus period preceding the TMS pulse (see Time and frequency analysis).The amplitude fluctuations of the components were visualized as an envelope of the filtered signals (figures 4(C) and (D)).Although the amplitude of the components varied throughout the 1.25 s period before the pulse, there were no consistent changes in the signal immediately before the stimulation onset.Notably, the amplitude of the component was higher in the trials from the congruent condition (i.e. the 'High' component in the 'High' trials) as opposed to the incongruent one (i.e. the 'High' component in the 'Low' trials), and this distinction was more pronounced with the 'High' component (figures 4(C) and (D)).To further investigate the significance of the signal's proximity to the TMS onset, we repeated the analysis using different time windows relative to the TMS pulse (figure 6).We observed a marginal gradual increase in accuracy with increasing proximity to the stimulation onset (figure 6(A)), from 66% ± 8% in the earliest window (1.25-0.75s) to 68% ± 8% in the latest window (0.5-0 s).The scalp topographies remained similar across the different windows (figure 6(B)).Within the 'High' patterns, there was a gradual shift in amplitude 'bridging' the left central-parietal, left frontal-central and right frontal regions.Similar to the analysis conducted on different frequency bands described earlier, the spatial patterns were averaged across all subjects included in the analysis.

Replication of the results
To validate our findings, we applied the same analysis to a different EEG-TMS dataset (see Validation analysis).The average classification accuracy across 11 analyzed subjects was 65% ± 5%, with 10 out of 11 subjects reaching statistical significance in classification accuracy with the threshold set at 59%.This accuracy is slightly lower but comparable to the main results (68% ± 8%).The group-average topography of the pattern was also consistent with the patterns of main analysis for both the 'High' and 'Low' conditions (figure 7).The 3% reduction in classification accuracy compared to the main analysis may be due the analysis algorithm being overfit to the main dataset, due to a lower number of electrodes in the EEG layout, or due to other possible differences in data collection.Nevertheless, these results successfully replicated our initial findings.

Patterns of spontaneous cortical oscillatory activity predict corticospinal excitability
We showed that aCSP components, derived from spontaneous oscillatory activity in the pre-stimulus  EEG signals, can predict the post-stimulus MEP amplitude.The variance of the aCSP components, reflecting their power in each pre-stimulus window, was predictive of the stimulation outcome in 95% of the analyzed subjects, with an average prediction accuracy of 68%.The achieved prediction accuracy is comparable to other ML-based approaches in decoding corticospinal excitability states from EEG-TMS data [15,16].This suggests that aCSP, guided by the readout of stimulation, revealed patterns of cortical activity that are relevant to the state of corticospinal excitability and, as such, may modulate the effect of TMS.
While we isolated source activity corresponding to high and low excitability states separately, the high-excitability component was particularly wellisolated (figures 4(A) and (C)).Compared to the low-excitability component (figures 4(B) and (D)), it was more distinct between the two experimental outcomes.This observation indicates that high MEP amplitudes have a more explicit relationship with the state of cortical activity, while low MEP amplitudes emerge for various reasons and may not be so clearly predicted by any singular oscillatory process.It is important to note here that high-and low-excitability components were used together in classification, so the classification success does not reflect their individual predictive power.
The achieved classification accuracy is relatively modest compared to the success of CSP application in other domains, such as brain-computer interfaces [21].There are several possible reasons for the mis-prediction of single-trial MEPs.In particular, MEP amplitude reflects excitability not only of corticospinal neurons in motor cortex but also of motoneurons in spinal cord [3,6,22].This leads to two implications: on one hand, corticospinal excitability as proxied by MEP amplitude may not be decoded from EEG signals in a deterministic way, and on the other hand, MEP is a noisy measure of the state of corticospinal excitability, susceptible to errors in labelling.Additionally, spontaneous oscillatory activity in EEG typically exhibits low signalto-noise ratio on short timescales, making it challenging to separate weak signals from background noise, even with advanced signal separation techniques [23].While these confounds are to some extent inherent to the classification of TMS-probed spontaneous states, they can still in principle be avoided by targeting other cortical states which are (1) clearly separable and (2) have a straightforward relation to the readout measure.Task-related states rather than spontaneous states could be one such target.

Individual scalp topographies are physiologically valid but not ubiquitous
The aCSP components were derived and selected individually, with the intention of isolating spatial patterns specific to each subject.Even if the underlying source activity was shared among individuals, differences in head geometry and EEG electrode placement would result in variations in the scalp distribution of the spatial filters and patterns.Despite the expected variability, we aimed to assess overlaps in the scalp topographies of the predictive components.If the most predictive spatial patterns exhibit similarity across individuals, it would indicate physiological validity as well as generalizability of the underlying neurophysiological phenomenon at the population level.
The distribution of the pattern magnitude across the scalp was significantly correlated across individuals for both high-and low-excitability patterns (see Spatial patterns).However, although the selected patterns exhibited greater similarity compared to a set of randomly generated topographies, their similarity was not higher than that of a random selection of other less predictive patterns generated by the aCSP on the same data.This suggests that, while the aCSP patterns were anatomically meaningful, there was insufficient evidence to conclude that the selected patterns shared a common source.This may be attributed to our selection of spatial patterns for further interpretation based on their eigenvalues, which represent the achieved level of separability between the two conditions [24,25].The same spatial pattern may be present in some or all of the subjects but have a relatively small eigenvalue for some of them, if patterns with stronger separability are available.Thus, we do not conclude that the predictive patterns vary across individuals, but rather that the most separative ones do.It is important to note that, while we selected only single components with the largest eigenvalues for further exploration of their spatial, spectral, and temporal features, several components from each experimental condition were sometimes combined in the classification to achieve the best prediction.

Spatial, spectral and temporal characteristics of the predictive components
The aCSP components obtained from EEG signals can be viewed as reconstructed oscillatory source activity, allowing for characterization in terms of their spatial localization, spectral composition, and temporal dynamics.
The high-excitability component was primarily localized in the left central-parietal scalp region, posterior to the stimulation site of the left motor cortex, and in the left frontal-central, right frontal and occipital areas (figure 3(A)).The alpha-specific oscillations dominated the component's spectral composition (figure 4(A)) and originated primarily in the left central-parietal region, posterior to the stimulated motor cortex (figure 5(B)).These observations suggest that the high-excitability component partially represents the sensorimotor mu-rhythm in M1 or the primary somatosensory cortex (S1).Indeed, phase and power dynamics of the sensorimotor mu-rhythm have been associated with changes in corticospinal excitability, although with inconsistent results [26][27][28][29][30][31][32][33][34][35][36].The beta-frequency range activity was less prominent in the component's spectral composition compared to alpha-oscillations (figure 4(A)), and its localization was distributed between left frontal-central, right frontal and occipital areas (figure 5(B)).Given the localization in the left frontal-central region anterior to that of the alpha-specific topography, this component may partially represent the sensorimotor beta-rhythm, propagating along anterior-posterior axis [19,37].Previous studies have associated the phase or power dynamics of the sensorimotor betarhythm with changes in MEP amplitude, although with varying findings [27,36].It may also simply be a harmonic of the mu-rhythm [38].Indeed, magnitude distribution of the broadband spatial pattern in our study (figure 3(A)) resembled a superposition of the alpha-and beta-specific patterns (figure 5(B)).
The phase shift progression observed in individuals sharing this pattern (figure 3(E)) suggests a travelling wave or phase-coupling along an anteriorposterior direction, possibly between M1 and premotor cortex or S1 and M1 [39][40][41].This observation is supported by the evolution of magnitude distribution of the high-excitability patterns with temporal proximity to the pulse (figure 6(B)).The observed shift in localization may be a temporal manifestation of a travelling wave between frontal and central-parietal regions.Despite these indications of the phase-shifted activity presence in the aCSP patterns, we did not verify its importance for the MEP prediction.The real CSP analysis, limited to only non-phase-lagged amplitude dynamics, predicted the MEP class equally well.Indeed, the variance of the real CSP and aCSP components in the same trials was significantly correlated.This indicates that taking phase-shifted activity into account did not provide any essential information for MEP amplitude prediction.Thus, another physiological interpretation of the observed high-excitability pattern could be a change in the dipole orientation rather than location, representing the spread of activation within M1 or S1 [42].
The low-excitability component was localized in the right parietal-occipital, medial frontal and bilateral temporal areas (figure 3(B)).The alpharhythmic activity prevailed in the power spectrum (figure 4(B)) and was localized in the medial parietal-occipital and frontal areas (figure 5(B)).The parietal-occipital alpha-rhythm is generally associated with closed-eyes state and idle state [43].Strigaro et al [44] found no effect of the eyes-open versus eyes-closed condition on either RMT or effective connectivity between visual and motor cortex.Moreover, our participants were instructed to fixate the eyes on a cross during the experiment.Still, there has been evidence of a modulatory connection between visual and motor cortex, probed with either visual or magnetic stimulation of visual cortex [44,45].Alternatively, the appearance of the effect of occipital alpha activity on MEP amplitudes could be explained by the coexistence of occipital alpha-and sensorimotor mu-rhythms within the same frequency range.The occipital alpha signals may be detectable during instances of low sensorimotor mu-activity, thus corresponding to lower MEP amplitudes, and vice versa (conceptually similar to the hand vs. foot imagery scenario in Blankertz et al [14].The phase shift progression suggests a possible underlying travelling wave or phase-coupling between parietaloccipital and frontal regions (figure 3(F)).However, this observation was not supported by further analysis (invariability of the spatial pattern across different time windows (figure 6(B), comparable prediction success of real and aCSP, and significant correlation of variance of real CSP and aCSP components across trials).
The observed spectral shape of the components suggests that alpha-frequency band activity may be relevant for both excitability conditions.Alpharhythm has been associated with top-down mechanisms of selective inhibition and information-gating [43].In contrast, the sensorimotor mu-rhythm has been suggested as a mechanism of temporally constrained facilitation, rather than inhibition [32].However, it is also possible that the alpha-peak reflects the dominant frequency band within the unfiltered sensor signals, either due to properties of the underlying sources or due to the general 1/f shape of the EEG spectrum [46].The latter is supported by our observation that prediction success depended on the spectral width of the bandpass-filter rather than on the specific choice of retained frequencies (figure 5(A)).Still, the sources of the high-and low-excitability components were clearly distinct, since their topographies were localized differently (figures 3(A) and (B)).
There were no systematic amplitude dynamics in either components time course within 1.25 s before the stimulation (figures 4(C) and (D)) and EEG signals from any latency within that period were equally successful in MEP prediction (figure 6(A)).Thus, a longer temporal window of the signal is necessary to verify the timescale of the neural activity involved.For reference, Hussain and Quentin [16] also found no difference in the success of MEP prediction when applying LDA to the power of oscillatory signals in different time windows within a 3-s period before TMS.

Other approaches to decoding corticospinal excitability with ML
It is worth briefly considering the difference between the current study and previous studies that used ML to decode corticospinal excitability from pre-stimulus EEG.Metsomaa et al [15] employed data-driven individual spatial and temporal filtering of EEG signal for decoding MEP amplitude from pre-stimulus oscillatory source activity phase-locked to the stimulation onset.The essential difference of our approach is that the targeted activity is not time-(or phase-)locked to the stimulation event.Similar to that study, the component does not necessarily originate from an anatomically restricted neuronal generator; instead, it may represent functionally coherent activity of a distributed network [13,25].However, due to the involvement of phase lags, the component's time course may aggregate not only simultaneous but also phase-delayed activity of the network.Hussain and Quentin [16] employed power of different spectral bands in the pre-stimulus EEG signals as predictors of the MEP amplitude in LDA.In the main analysis of the current study, the EEG signals were broadband-filtered and thus were not frequencyspecific.Moreover, rather than using sensor-level power features, we reduced signal dimensionality with aCSP and then supplied the isolated components to the classifier.

Limitations of the study
The use of MEP amplitudes as means for categorizing stimulation outcomes into discrete conditions has its specifics.MEP amplitude, like corticospinal excitability, is inherently a continuous measure rather than a discrete one [47].Discretizing the continuous measure, while necessary for the current analysis, omits a certain level of the underlying complexity of the data.Nonetheless, we opted for the discrete measure for two main reasons.Firstly, in the simplest brain statedependent stimulation paradigm, the decision space for stimulation is discrete (deliver or not deliver).Thus, having discrete information about the ongoing brain state (e.g.present or absent) makes it easier to inform the decision.Secondly, when applying statistical analyses to continuous measures (e.g.regression), assumptions need to be made about the shape of the relationship between the analyzed variables (e.g.linear or non-linear).Nevertheless, future studies may consider incorporating non-discrete readout or extracting non-discrete brain states using other signal extraction methods besides CSP.
When interpreting the results of our analysis performed on the signals filtered in different frequency ranges, it is important to consider that the use of narrowband frequency filters may distort the signal, worsening the estimation of spatial covariance [48].Indeed, we observed lower classification accuracy when narrower-filtered signals were analyzed (see Time and frequency analysis).However, adopting a more frequency-specific approach seemed to improve the isolation of the aCSP components in terms of their spatial localization (figure 5(B)).Future studies may go for either one of the two approaches depending on whether the focus is on maximization of prediction accuracy or maximization of interpretability of the results.

Conclusion
We employed a machine-learning approach in combination with blind-source separation in the form of aCSP to derive predictors of corticospinal excitability from spontaneous EEG activity.The isolated oscillatory patterns represented network-level oscillatory activity, and the variance of these patterns predicted the MEP amplitude.We found predictive activity within the analyzed 0.5 s time window before the TMS pulse in the 8-30 Hz frequency range.The activity predictive of high corticospinal excitability was localized in the lateral central-parietal region close to the stimulated motor cortex.The activity predictive of low corticospinal excitability was localized in the medial parietal-occipital and frontal areas.The predictive components from both conditions had a spectral peak in the alpha-frequency band.Overall, we established a data-driven approach to uncovering network-level oscillatory activity that modulates TMS effects.The aCSP approach requires no anatomical priors, while being physiologically interpretable, and can be employed in both exploratory investigation and brain state-dependent stimulation.

Figure 1 .
Figure 1.Analysis procedure.(A) For each subject, all experimental trials were separated into high and low excitability conditions according to their corresponding MEP amplitude.(B) Preprocessed pre-stimulus EEG signals from 0.5 s before each TMS pulse were then divided into two groups according to the condition labels.(C) Covariance matrices of the EEG signals from both conditions were averaged across trials for each condition and used in aCSP analysis.(D).The aCSP produced spatial filters aimed at isolating signal components that maximize separation between the two experimental conditions.The effectiveness of separation was tested in the following way.(E).The pre-stimulus EEG signals were spatially filtered with the aCSP filters.(F).The variance of the filtered signal components in each trial was employed as features in LDA classification to predict the excitability condition.The classification accuracy was measured as a proportion of correctly classified trials with respect to the original labelling based on the MEP amplitude.With the exception of the condition labelling stages (A)-(B), the aCSP analysis underwent a 5-time 5-fold CV procedure (outer CV layer, (C)-(F)).The average classification accuracy across all CV folds of this layer was taken as an overall classification accuracy of the given subject.The hyperparameters used in aCSP were derived via an additional 5-fold CV procedure on each iteration of the outer CV layer (inner CV layer, (C)-(F)).

Figure 2 .
Figure 2. Example results from a single subject.(A) Separation of trials into the 'High' (in red) and 'Low' (in blue) corticospinal excitability conditions based on the MEP amplitude.Black line corresponds to a moving median of 150th order separating the data into the two conditions.(B) Spatial filters generated with aCSP for the 'High' (top) and 'Low' (bottom) conditions.Complex-valued spatial filters are visualized as pairs of magnitude (absolute values, left) and phase (right) topographies.(C) Spatial patterns derived from the spatial filters (in (B)).The complex-valued spatial patterns are visualized as pairs of magnitude (left) and phase (right) topographies.Phase topographies depict spatial distribution of phase lags across the scalp with respect to a reference channel (indicated with a black dot).(D) Dynamical depiction of spatial patterns for the 'High' (middle row) and 'Low' (bottom row) conditions as a function of phase (top row).The spatial patterns from (C) can be alternatively visualized as the change in voltage distribution across the scalp as a function of a phase of an oscillatory signal projected onto the scalp.See also figure S1. (E) Time course of the 'High' (top) and 'Low' (bottom) aCSP components.Time courses are visualized as the real part (left) and the envelope (right) of the filtered signal, averaged across trials.For comparison, both components were extracted from the EEG signals from both 'High' (in red) and 'Low' (in blue) trials.Dashed vertical lines indicate the time of the TMS pulse.F. Power spectrum of the 'High' (top) and 'Low' (bottom) aCSP components.Power spectrum was calculated on the spatially filtered signals.

Figure 3 .
Figure 3. Spatial patterns.(A) and (B).Group-average magnitude topography for the 'High' (A) and 'Low' (B) spatial patterns.The 'High' pattern corresponds to a component explaining the most variance in the high-excitability trials and least variance in the low-excitability trials, while the 'Low' pattern corresponds to a component explaining the most variance in the low-excitability trials and least-in the high-excitability trials.The colors represent the distribution of amplitude of the component across the scalp.The average topographies were calculated on the subjects with significant prediction accuracy (N = 19).See also figures S2 and S3.(C) Group-average magnitude topography for the 'High' spatial pattern calculated on a subset of subjects with a similar 'High' topography (N = 6).(D) Group-average magnitude topography for the 'Low' pattern calculated on a subset of subjects with a similar 'Low' topography (N = 6).(E) Group-average phase topography for the 'High' pattern calculated on a subset of subjects with a similar 'High' magnitude topography (N = 6, same subset as in (C)).The colors represent the phase shift (in degrees) with respect to a reference channel (FCC3h channel, indicated with a black dot).All phases were subtracted from the phase in the reference channel.The channels exhibiting the magnitude below the median of the magnitude distribution were masked (in grey).(F) Group-average phase topography for the 'Low' pattern calculated on a subset of subjects with a similar 'Low' magnitude topography (N = 6, same subset as in (D)).The colors represent the phase shift (in degrees) with respect to the reference channel (Cz channel).

Figure 4 .
Figure 4. Power spectrum and amplitude time course of the aCSP components.Spatial filters were computed on the signals within a 0.5 s window before the TMS pulse (dashed vertical line in (C) and (D)) and applied to a longer 1.25 s window for visualization.The filters were applied to a set of trials that were not used in their computation.'High' and 'Low' labels represent which excitability condition either the EEG signals or the components correspond to.(A) and (B).Power spectrum of the 'High' (A) and 'Low' (B) components isolated from the EEG signals recorded during the 'High' (in red) and 'Low' (in blue) trials (mean ± standard error of the mean (SEM) across subjects).Power spectrum was averaged across trials and then across subjects with significant classification accuracy (N = 19).(C) and (D) Time courses of the amplitude are visualized as the envelope of the 'High' (C) and 'Low' (D) components isolated from the signals from the 'High' (in red) and 'Low' (in blue) trials (mean ± SEM across subjects).The time courses were averaged across trials, normalized to unit-norm for visualization and then averaged across subjects with significant classification accuracy (N = 19).

Figure 5 .
Figure 5. ACSP analysis in different frequency bands.(A).Classification accuracy of aCSP analysis of EEG data filtered in four different frequency bands (mean accuracy ± SEM across subjects, N = 20).(B).Magnitude topographies of the 'High' (top) and 'Low' (bottom) spatial patterns in the four analyzed frequency bands.The 'High' and 'Low' labels represent, which excitability condition the patterns correspond to.The colors represent the distribution of the component's amplitude across the scalp.The average topographies were calculated on all analyzed subjects (N = 20).Of note, higher accuracy in the beta-and low gamma-band is at least partially due to broader bandpass filtering of the EEG signals, resulting in less signal distortion.

Figure 6 .
Figure 6.ACSP analysis in different time windows.(A).Classification accuracy of aCSP analysis of EEG data from different overlapping time windows in the pre-stimulus period (mean accuracy ± SEM across subjects, N = 20).(B).Magnitude topographies of the 'High' (top) and 'Low' (bottom) spatial patterns in the four analyzed time windows.The 'High' and 'Low' labels represent, which excitability condition the patterns correspond to.The colors represent the distribution of the component's amplitude across the scalp.The average topographies were calculated on all analyzed subjects (N = 20).

Figure 7 .
Figure 7. Spatial magnitude patterns from the validation dataset.Group-average magnitude topographies for the 'High' (left) and 'Low' (right) components.'High' and 'Low' labels represent which excitability condition the patterns correspond to.The colors represent the distribution of the component's amplitude across the scalp.The average topographies were calculated on the subjects with significant classification accuracy (N = 10).