MEP and TEP features variability: is it just the brain-state?

Objective. The literature investigating the effects of alpha oscillations on corticospinal excitability is divergent. We believe inconsistency in the findings may arise, among others, from the electroencephalography (EEG) processing for brain-state determination. Here, we provide further insights in the effects of the brain-state on cortical and corticospinal excitability and quantify the impact of different EEG processing. Approach. Corticospinal excitability was measured using motor evoked potential (MEP) peak-to-peak amplitudes elicited with transcranial magnetic stimulation (TMS); cortical responses were studied through TMS-evoked potentials’ TEPs features. A TMS-EEG-electromyography (EMG) dataset of 18 young healthy subjects who received 180 single-pulse (SP) and 180 paired pulses (PP) to determine short-intracortical inhibition (SICI) was investigated. To study the effect of different EEG processing, we compared the brain-state estimation deriving from three published methods. The influence of presence of neural oscillations was also investigated. To evaluate the effect of the brain-state on MEP and TEP features variability, we defined the brain-state based on specific EEG phase and power combinations, only in trials where neural oscillations were present. The relationship between TEPs and MEPs was further evaluated. Main results. The presence of neural oscillations resulted in more consistent results regardless of the EEG processing approach. Nonetheless, the latter still critically affected the outcomes, making conclusive claims complex. With our approach, the MEP amplitude was positively modulated by the alpha power and phase, with stronger responses during the trough phase and high power. Power and phase also affected TEP features. Importantly, similar effects were observed in both TMS conditions. Significance. These findings support the view that the brain state of alpha oscillations is associated with the variability observed in cortical and corticospinal responses to TMS, with a tight correlation between the two. The results further highlight the importance of closed-loop stimulation approaches while underlining that care is needed in designing experiments and choosing the analytical approaches, which should be based on knowledge from offline studies to control for the heterogeneity originating from different EEG processing strategies.


Introduction
The use of non-invasive brain stimulation has expanded in the last years, leading to new technical developments, protocols, and therapies.Along with innovations and improvements, the wide open-loop use of these techniques has also raised further questions on their working principles as a high heterogeneity in results and small-to-medium effect sizes have been highlighted [1][2][3].The ongoing brain oscillatory activity has been suggested as one of the main causes of such variability [2,4,5].It is just in recent years that researchers have begun to study behavioural and electrophysiological responses to stimulation in relation to ongoing cortical activity (table 1, offline protocols, top rows); and as more studies hinted to an effect of oscillatory activity, the first brainstate-dependent stimulation protocols started to be published (see table 1, online protocols, bottom rows).
The most-widely used techniques for such experiments are transcranial magnetic stimulation (TMS) in combination with electroencephalography (EEG), and electromyography (EMG).The first allows to perturb endogenous oscillations.The second permits to both read and predict the brain state and record the cortical response to the stimulation (i.e.TMS-evoked potentials-TEPs).The third allows to evaluate corticospinal excitability by studying motor evoked potentials (MEPs).Figure 1 depicts the three technologies and how they may be used together to perform and evaluate brain-state-dependent cortical stimulation.
In previous studies, the brain state has been defined through the estimation of power spectral density (PSD) and phase of ongoing oscillations (from now referred to just as 'phase') of the sensorimotor system covering the alpha (or mu, 8-12 Hz) and the beta (13-30 Hz) bands.Specifically, the PSD relates to the size of the neural population that is firing synchronously [7,8], whereas the phase rather represents the state of the population (e.g.excitation or inhibition) [9][10][11].
The introduction of the brain state to explain MEP variability did not convey much homogenization-see results from table 1.For example, when looking at the relationship between MEP amplitude and PSD in the alpha band when subjects are at rest, 3 out of 14 offline studies (i.e.studies where trials are sorted retrospectively according to the brain-state at the stimulation) reported a negative relationship, 7 of those did not find any effect and the remaining ones found a positive relationship.When extending to online studies (i.e.studies where the stimulation was triggered according to the ongoing brain-state as measured with EEG), 2 out of 7 experiments found a negative relationship, 4 a positive one, and the remaining ones no effect.Only very few online studies have been performed (N = 10) and of these, 60% came from the same research group, where the same inclusion criteria and processing methods have been applied-importantly, reproducibility within this group was high.The discordant outcomes between groups were often attributed to stimulation protocol differences, such as number of trials and subjects, stimulation intensity, and interstimulus interval (ISI).In addition to that, all the signal processing steps, involving both data cleaning and filtering, have been shown to play a critical role in phase accuracy estimation [12], but have not yet been studied in terms of studies' outcomes, such as MEP variability.As outlined in table 1, studies use different lengths of pre-stimulus window, spatial and frequency filters, and power estimation methods, and have often not reported the exact processing, especially in earlier work.Finally, the statistical methods used for comparison also differed across studies.
Significantly fewer experiments have focused on the relationship between brain state and cortical responses through TEPs (see table 2).This can be explained by the multidimensionality of TMS-EEG induced signals and the rather recent conceptualisation of this technique.As for the MEP results, heterogeneity is still present.Importantly, unlike MEPs where the amplitude and latency of one deflection are looked at, in the case of TEPs, multiple peaks can be investigated as well as other features such as the mean field potential, inter-trial coherence (ITC), and evoked power, to name a few.In the literature, different studies have focused on different features but the most-present features are related to the amplitude of the well-known TEPs components at the group level (i.e.P25, N45, P60, N100).For instance, N100, a well-accepted marker of intracortical excitatoryinhibitory network activity, seems to be modulated by the phase of the brain state across various works [13][14][15]; yet, N100 has been recently suggested not to be related to the motor network excitability [16].
The present work results as a consequent followup of our previous analyses focused on the effects of EEG processing on phase prediction accuracy [12].Here, we provide first and strong evidence of the impact of EEG processing on the outcome of brainstate-dependent stimulation as determined by MEP and TEP features.At this regard, the effect of presence and absence of neural oscillations on the preprocessing will be considered as well, as it has recently been suggested to be of critical importance in EEG spectral analyses [6,17].Moreover, the present work contributes to the current state of the art regarding the understanding of motor cortex physiology in the alpha band.

Data acquisition
The dataset was acquired in the framework of the TiMeS project, of which the full protocol can be found at [47].TMS-EEG experiments were run with 18 healthy young right-handed participants (7 females, 27.0 ± 2.8 years old).All participants signed an informed consent in accordance with the Declaration of Helsinki and a sheet related to the risks of TMS.The study was approved by the ethical committee of the Canton of Vaud, Switzerland (no.2018-01355).Other aspects of these data have been published in [48,49].
The experiment was performed in a Faraday cage, the participants were seated in a comfortable chair, with elbows flexed at 90 • and hands relaxing on a table ; with eyes open, participants were instructed to Figure 1.A set-up for brain-state-dependent cortical stimulation.For brain-state-dependent cortical stimulation, the subject is sitting comfortably on chair at rest.The brain activity is recorded with EEG and is processed by a computer (box on the right); among the main processing steps there is spatial filtering, the calculation of PSD, here also defined as the presence or absence of neural oscillations (as suggested by [6]), frequency filtering and an algorithm to forecast the phase of interest.Whenever this time point is reached, cortical stimulation by TMS is triggered and electrophysiological responses may be studied at the muscle level with EMG (here placed on a hand muscle) or at the cortical level (see figure 2).look at a fixation cross in front of them.The EEG system consisted of 64 Ag/AgCl TMS-compatible electrodes in a 10-20 system (EEG BrainCap-MR BrainVision LLC, North Carolina, USA) with 5 kHz sampling rate and 1 kHz high cut-off.Synchronously, muscle activity was recorded from the left first dorsal intraosseous.EMG was captured with a pair of disposable Ag-AgCl electrodes; the signal was amplified and sampled at 3 kHz using a Noraxon DTS Receiver (Scottasdale, Arizona, United States) with the bandpass filter at 10 Hz to 1000 Hz (analog Sallen-Key for high-pass filter and digital FIR filter with order 128 for the low-pass; the gain was set at 500), and digitized at 5 kHz using Signal software (Cambridge Electronic Design Limited, Cambridge, UK).
TMS was delivered with a MagPro X100 stimulator connected to an MC-B70 coil (Magventure, Farum, Denmark) at rest.First, hot-spot hunting was conducted manually, following standard guidelines (as originally proposed in works such as Rossini et al [50] and Van de Ruit et al [51].This process involved exploring the cortical motor map of the targeted muscle, starting from an initial point over the motor hand knob, and searching for the target that elicited the largest and most stable MEPs.The hotspot was then defined as the coil location and orientation delivering the highest MEP peak-to-peak amplitude.Along the experiment, a neuronavigation system (Localite neuronavigation software, Localite GmbH, Germany) was exploited to register the hotspot and all the coil positions at stimulation.The coil was rotated at 45 degrees with respect to the midline and stimulation was delivered with a posterior to anterior current.
Participants underwent six stimulation blocks with 60 pulses each, with two alternating pseudorandomised sequences containing monophasic SP and PP using SICI conditions (n = 30 for one block, 180 in total, for each condition).The latter were programmed as two pulses with 3 ms inter-pulse interval with the conditioning stimulus at 80% of the resting motor threshold (rMT) and the test stimulus at test intensity (TI).RMT was defined as the minimal intensity to evoke 5 out of 10 MEP with peak-topeak amplitude larger than 0.05 mV [50]; TI was defined as the minimal intensity giving consistently a peak-to-peak amplitude of at least 0.5 mV, if this was not found, 120% of the rMT was used.The interstimuli ISI was of 4 s with a 25% jitter.During the stimulation blocks, participants were asked to stay relaxed and keep their eyes open and look at the fixation cross.Participants were wearing noise-cancelling headphones playing white noise to cover the TMS sound and a thin layer of foam was applied on the coil to reduce the acoustic interference.The full length of the experiment was of about 2.5 h; see supplementary material for further information on the overall experimental protocol.

EEG data analysis
Visual inspection was performed to manually remove and interpolate bad trials and channels.EEG processing was done using EEGLAB [52], TESA [44], and Fieldtrip [53] in Matlab 2021.b (Mathworks Inc., Table 1.Synoptic table of literature review regarding articles where MEP features (amplitude and latency) have been considered in relation to the brain state at the time of TMS stimulation.Both offline and online studies are reported; the type of protocol is highlighted in bold in the column 'Experiment condition' .In this table we report in a summary the following information: reference paper and subjects demographics (i.e.sample size, age, and additional information); experiment condition (i.e.protocol with information of activity performed by subject, TMS related parameters and EEG and EMG electrodes used); EEG processing pipelines including both preprocessing, choices of hyperparameters (e.g.window duration of epoch) and how PSD and phase have been estimated); results in the alpha and beta bands to related MEP features in terms of phase (PH), PSD or interaction effect (PH x PSD); statistics used.If an analysis was not performed by the article, results are reported as NA (not applicable); if a method was not specified, it is reported as NR (not reported).Abbreviations and units used going from left to right and top-bottom: N is the sample size, F is the number of female subjects; age is reported in years, either mean  Maki & Ilmoniemi, 2010 [22] healthy right-handed,   Table 1.

EEG preprocessing and brain-state definition
To estimate the brain-state, we employed the following approach: the pre-cleaned data was epoched into trials starting 0.75 s prior to the TMS pulse until the TMS pulse [12].The trials were spatially filtered at the C3 electrode level with a Hjort-Laplacian (i.e.remove the mean of the neighbours FC1, CP1, FC5, CP5) [55].This filter is widely applied when studying motor cortex signal at the electrode level and we believe it is a good approximation for all subjects for defining the left hand-knob area activity, although not necessarily always the closest electrode to the stimulation hotspot.Power was estimated using the Welch approach and by summing the magnitude of the spectra in the whole frequency band of interest.
From the power spectra, we then evaluated if a peak was present in the band of interest using the approach and code proposed by [6] with default parameters.Thus, a peak in the frequency band of interest (i.e.8-12 Hz) was reported if a peak was actually present in the periodic component of the spectra and if its value was larger than the aperiodic component.The binary output was used to define whether a neural oscillation was present in each trial.Finally, phase extraction was done following the processing found previously by a grid-search [12]: epoch length of 0.75 s prior to the TMS pulse, 2nd order bandpass Bessel filter with a bandwidth of 1 Hz around main frequency of the epoch.The filter was applied after zero-padding the epoch of the same length of the original signal; after filtering, 30% of the end of the signal was cut and the phase at the stimulation point estimated by using as period the double the distance between the last peak and trough found and as starting phase the phase found with the Hilbert function [56].In addition to the pipeline described above, we have reproduced the methodology of [30,35,39] and [41] to test their results on our dataset and evaluate the importance of EEG processing choices on outcome measures.According to the estimated phase in radians between −π and +π, this was categorized into four groups following a sine wave: trough (−3π/4-−π/4), rising (−π/4-π/4), peak (π/4-3π/4) and falling (>3π/4 or <−3π/4).PSD was kept as a continuous numerical variable when studying the effects on MEP features, whereas we categorized it into quartiles per subject and time point to have the very low (below 25th percentile), low (between 25th-50th percentile), high (between 50th-75th percentile), and very high (above 75th percentile) classes to study TEPs.The difference in approaches lies in the difference between MEPs (i.e.single-trial response) and TEP (i.e.averaged response).Specifically, PSD was not categorized for MEPs models to avoid having too few trials per interacting groups.

TMS-evoked potential features
TMS-EEG data were analyzed with MATLAB (The MathWorks, USA) and preprocessed using the EEGLAB [43] and TESA [45] toolboxes, following the pipeline suggested by [44] and reported in [48].First, raw EEG data were epoched and demeaned between −0.5-1 s around the TMS pulse, the signal between −5 and 20 ms being removed and interpolated (using cubic interpolation).The signal was then downsampled to 1 kHz, and the previouslyinterpolated period replaced by constant amplitude data.Bad channels and epochs were visually identified and removed.A first round of independent component analysis (ICA) was run to remove large-amplitude artifacts, using FastICA and the standard parameters proposed in TESA.Data around TMS pulse were re-interpolated before signal filtering between 1 and 100 Hz.Then, the previously-interpolated period was replaced by constant amplitude data before running a second run of ICA to remove all remaining components linked to pulse, muscle, decay and ocular artifacts from the data.Finally, missing channels and constant amplitude data around TMS pulse were interpolated, and the clean EEG signal was re-referenced using average reference.TEPs were computed by averaging clean trials for each subject.Once TEPs were obtained, the single trials were labelled according to the feature of interest (i.e.categorial phase, categorical PSD, or MEP size), averaged per subject and time point.Numerous features from time and frequency domain were extracted, summarized in table 3 and in figure 2.

EMG data analysis
EMG trials were rejected if the standard deviation of the baseline (400 ms before TMS onset) was bigger than 0.05 mV.MEP peak-to-peak amplitude was computed as the difference between the maximum and minimum of the EMG signal in the window 10-50 ms after the stimulation.The latency was computed with the method published in [59].

Statistics
Effects of presence of neural oscillation, phase, and power of the brain state in the alpha band were studied using generalized linear mixed effect models (glmer function in R [60]) with independent fixed variables of phase (categorical), power (continuous) and their interaction; and random intercept given by the subject.The dependent variable was the peak-to-peak amplitude or latency of the MEPs.Only trials with neural oscillation and where neither EMG nor EEG artefacts were present were considered.Different analyses were run for the two pulse conditions.Post-hoc analyses used estimated averaged means (emmeans package in R).The supplementary information provides important details on how specific GLMM were chosen, together with specific transformation of variables and how results can be  3.
interpreted.A GLMM with gamma distribution was used to describe the MEP amplitude and a log-linear model to describe the MEP latency.
Investigation of the brain-state effect on TEP features reported in table 3 was also performed using linear mixed effects models, with post-hoc estimated marginal means with Bonferroni correction for multiple comparisons.Pulse conditions were evaluated with separate models.

Dataset
On average there were 178.3 ± 7.0 trials per subject and stimulation type (both for PP and SP); 86% of trials were retained after EMG and EEG artefact rejection.The EMG preprocessing removed 3% (5% for PP) of trials, mostly in one subject.An additional source of trials removal came from the binarization of trials according to the presence of neural oscillations.In the cleaned data, in the alpha band around 66% of trials showed detectable neural oscillation (i.e.presence of a peak in the frequency band of interest in the periodic part of the power spectra [6]), when the PSD was computed with the Welch method over the 0.75 s window.

Effects of EEG preprocessing on PSD and phase estimation
A first step aimed at analysing the effect of different EEG processing steps on phase and PSD estimation.We looked at the steps of spatial filtering (here designed as a Hjort, common average reference CAR, or no filter applied), PSD estimation Table 3. Features from TEPs.On the left column the group of features according to the analysis, on the right the respective specific features extracted.All features, except those for the time domain, were measured in the time region of interest (TOI) = 20-80 ms after stimulation.Only for the area of the TEP curve we looked at the longer TOI = 20-300 ms.Some graphical depictions of such features can be retrieved from figure 3. Abbreviations: GMFP = global mean field potential, LMFP = local mean field potential, ITC = inter-trial coherence, RQS = regression quality score.

Specific features TEP in time domain
Amplitude and latency of ideal peaks (P25, N45, P60, N100, P 180, N250); complexity index; first found peak amplitude and latency.Mean of and area under the TEP curve.For these features the TEP was extracted from the average of motor cortex channels: C3, C5, C1, FC3, CP3.
Global/Local mean field power (G/L MFP) GMFP area (first, second half and all); LMFP area; LMFP complexity index.The GMFP area was computed across the whole brain, whereas LMFP was computed for the motor cortex channels: C3, C5, C1, FC3, CP3 Time-frequency analysis Amplitude and latency of max and min power; mean power; zero-crossing.Time-frequency analyses was computed with a Morlet wavelet.Baseline (500-2 ms before the stimulation) correction was applied, by dividing by the mean of the baseline.For these features the TEP was extracted from the average of motor cortex channels: C3, C5, C1, FC3, CP3.
Inter-trial coherence (run for theta, alpha, low beta, high beta, beta, and all frequencies) Mean and maximum ITC and latency of maximum ITC.Time-frequency analyses was computed with a Morlet wavelet.Baseline (500-2 ms before the stimulation) correction was applied, by dividing by the mean of the baseline.For these features the TEP was extracted from the average of motor cortex channels: C3, C5, C1, FC3, CP3.
Regression quality score (RQS) [57,58] Feature that studies the overall evoked dynamics by regressing a reference TEP onto single-trials possibly coming from another class.It is described by the resulting regression t-score (the higher, the more similar the two curves).We talk about paired RQS (pRQS) when single-trials and reference TEP belong to the same class; in the other scenario we talk about unpaired RQS (uRQS).The RQS is computed using TEPs and the clean trials from which they were computed.The goal of the score is to compare the (dis-) similarity between the (averaged) local TEP (computed from C3, C5, C1, FC3, CP3) and the single-trials, here in the early part of the response (from +20 to +80 ms after TMS pulse).The comparison is done by linearly regressing the two curves, with the following approach: s j (t) = β × x i (t) + ε (t) , with t ϵ [20,80] ms Where i and j are {'peak' , 'trough' , 'rising' , 'phase'}, or {'very-low power' , 'low-power' , 'high-power' , 'very-high power'}, or yet again {'very-low MEP' , 'low-MEP' , 'high-MEP' , 'very-high MEP'} conditions, depending on the analyses that is being performed.The t-statistic associated with the TEP factor xi resulting from the regression were saved for each trial, comparison (i.e.phase, power, MEP size), and subject.For group-analyses the scores were averaged across subjects and trials.If needed, further details can be found in Passera et al 2022 [58].
(Welch, Fast Fourier transform FFT, autoregressive model using Burg) and overall approach (methods proposed here, in [30], and in [41]).We discovered that spatial filter selection was rather critical, especially concerning phase estimation (supplementary figures 2(A) and (D)), whereas PSD choice, though creating some differences, gave overall similar results achieving a high correlation for both PSD estimation method and overall approach (supplementary figures 3(A) and (D)).In terms of phase estimation, differences were important (supplementary figure 3(G)).
We then performed the same analyses by grouping trials according to the presence/absence of neural oscillations in the alpha band.Higher correlations between different approaches were obtained for spatial filter, PSD, and overall methodology choices on both PSD and phase estimation when taking trials with alpha oscillations (supplementary figures 2(B) and (E); supplementary figures 3(B), (E) and (H)) compared to trials without a neural oscillation (supplementary figures 2(C) and (F); supplementary figures 3(C), (F) and (I)).This result has however to be tempered, given that the number of trials without alpha oscillations was considerably smaller compared to trials with.Indeed, as it can be observed from the numbers in Supplementary figures 2 and 3, very few trials resulted being consistently absent of a neural oscillation across approaches.Given the higher accordance of results using neural oscillations, for the next analyses we first looked at the effect of presence/absence of neural oscillations on the variable of interest and then performed analyses only with trials with oscillations.

MEPs and brain alpha oscillatory activity
Brain-state effects on the MEP peak-to-peak amplitude, a proxy for corticospinal excitability, was studied using the methodology proposed in [12].For single subject distribution of MEPs and RMTs, see supplementary figure 4. Separate analyses were performed for SP and PP conditions.No significant effect of the presence or absence of neural oscillations was observed regarding MEP amplitude (permutation test, p = 0.480).
When considering only trials with neural oscillations, we found a significant effect of phase and PSD on the MEP amplitude during SP.Post-hoc analyses showed that MEPs elicited during the peak phase of the alpha band were significantly smaller than those elicited during the rising (p = 0.045, estim = −0.12)and trough phases (p = 0.034, estim = −0.13).The trough phase elicited also stronger MEPs compared to those of the falling phase, though only at a trend level (p = 0.093, estim = −0.10).There was an overall positive correlation between MEP amplitudes and PSD (p = 0.030).The strongest relationship was observed in the trough and rising trials.The results are depicted in figure 3(a) and the full statistics reported in table S1.During PP stimulation, the significant effect of phase persisted with the peak trials eliciting the smallest MEPs compared to the three other phases (p < 0.05 and estim < 0.14 for all comparisons).Although the main effect of PSD disappeared, the interaction effect of Phase X PSD was present, with trough trials presenting a negative MEP-PSD correlation (p = 0.026), which was significantly different than the trend in the rising and peak trials.A similar difference was observed for the falling label.When MEP latency was used as dependent variable, a similar effect of phase observed for the amplitude was seen during SP, with peak trials having longer latencies compared to trough and rising ones (p = 0.029 and p = 0.033, respectively); falling trials elicited longer latencies than trough ones, as a trend (p = 0.072).Latencies and amplitudes had a negative correlation regardless of the brain state (r = −0.405,p < 0.001 when tested with logarithmic correlation tested with Pearson's test).No effects were observed in the PP ondition.

Effects of brain-state extraction method on MEP
Given the contradictory results of some online brainstate-dependent stimulation experiments and the differences introduced by the EEG processing pipelines, we used our dataset to reproduce and compare the results of [30] and [41], here referred to as the Tübingen and Copenhagen methods following their main city of provenance.To allow for a better comparison across studies we divided the alpha wave phase into peak, falling, trough, and rising and no threshold was applied on the PSD.Given different windows of interest and PSD estimation methods, only 1229 trials were common to all the approaches, and among those, around 40%-45% were labeled with different phases between each other, with only around 1%-2% of trials being labeled in anti-phase when doing a pair-wise comparison of the approaches.When using the Tübingen method (figure 3(c)), we observed a main positive effect of PSD (p = 0.010) and of phase, with trough trials being significantly bigger than peak and falling (p = 0.005, estim = 0.16 and p = 0.027, estim = 0.11, respectively).The rising phase also elicited larger MEPs than the peak phase, as a trend (p = 0.068, estim = 0.10).The interaction effect was also present, with trough trials presenting a positive MEP-PSD relationship, which was significantly different from that of the other labels (p < 0.01).When the Copenhagen method was used (figure 3(b)), the main effect of PSD was not present, but a positive MEP-PSD relationship remained as a trend for the trough trials (p = 0.062), which also elicited significantly bigger MEPs compared to the falling (p < 0.001, estim = 0.27) and the rising (p = 0.019, estim = 0.16) phases.For consistency of comparison, we used the same generalized linear effect model for the three approaches.Note that results belonging to the original statistical model are reported in tables S6 and S8 for [30] and [41], respectively.As mentioned in the previous section, the effect of the statistical model is again evident, with the Copenhagen method losing all the significant effects, when the original method is employed.Importantly, although the number of trials labeled differently is significant, the results of the mixed-effect models share the phase effect with trough trials eliciting stronger MEPs and falling trials smaller ones.

TEPs and brain state
TMS coupled with EEG allows to study cortical responses to perturbations.Similarly to corticospinal responses, we investigated if the brain-state could affect TEP features (figure 2, table 3).Again, after studying the effect of presence/absence of neural oscillations, only trials with an alpha peak were considered.Given the considerable amount of features studied (table 3), we here report only the significant results.Complete results can be found in supplementary tables 9 and 10.
In the time-domain, the presence of neural oscillations related to a moderately lower global mean field potential (GMFP) in both SP (p = 0.017) and PP (p = 0.001) conditions compared to the trials without neural oscillations.From the time-frequency analyses, we observed a larger mean and maximum ITC in trials with neural oscillations (p < 0.001 for both features and pulse conditions).Brain-state effect on MEP peak-to-peak amplitude in single-pulse condition according to different approaches.The figure shows the results from the generalized linear mixed models for the alpha band (see tables S1, S5 and S7).On the left the main effect of phase is shown and, on the right, the main effect of PSD (black curve) and its interaction with phase (colored curves).The MEP amplitude used for the models was kept raw, but for consistency and easier view it was z-score normalized per subject for the plot.The PSD used in the model was z-score normalized per subject.(A) Lausanne approach-full statistical results are reported in table S1; (B) Copenhagen approach [41]-full statistical results are reported in table S7; (C) Tübingen approach [30]-full statistical results are reported in table S5.Significant differences are reported.If the main effect of PSD was observed, it is reported, together with its trend at the top left corner of the phase-PSD interaction plot.The estimate of the model is also reported.Color legend: trough phase: green, rising phase: orange, peak phase: violet, falling phase: pink.Significance abbreviations: p < 0.1: '.' , p < 0.05: ' * ' , p < 0.01: ' * * ' , p < 0.001: ' * * * ' .Effects of oscillation phase were visible during SP trials on the local mean field potential (LMFP, figure 4(b)), with peak trials eliciting a smaller potential than trough ones (p = 0.042).The same significant difference was seen for the P60 amplitude (p = 0.037; figure 4(a)).Significant results of alpha phase were seen for the regression quality score (RQS) feature where the paired RQS (pRQS) of the trough trials was significantly higher than the unpaired RQS (uRQS) of peak and rising (p < 0.001) and the pRQS for the falling phase was higher than the uRQS with the peak phase (p = 0.040).RQS trends and the phase effect on P60 were present also in the PP condition (figure 4(d)).PSD divided into quartiles per subject (very low, low, high, and very high) had effects on time-frequency analysis features: in SP the maximum evoked power was significantly bigger in very high trials compared to the other groups (p ⩽ 0.004) and significantly lower in the very low trials (p ⩽ 0.020, figure 4(c)).The maximum value of ITC was bigger for high power trials than very low ones (p = 0.023).This trend was seen also in PP, where in addition the mean ITC was significantly bigger in the high trials compared to both very low and very high ones (p ⩽ 0.026).In this condition, the mean evoked power was significantly smaller for the very low power trials compared to the other groups (p ⩽ 0.002) and the mean TEP was bigger in very high compared to very low power trials (p = 0.053).

TEP and MEP correlation
MEP peak-to-peak amplitude and TEP features have been previously found to be correlated [16,58].We sorted TEP trials according to the MEP size (quartiles per subject) and studied if any relationships were present; analyses were performed separately for SP and PP conditions.We here report only significant results, complete results of statistics can be found in supplementary table 11.MEP size correlated with the P60 amplitude, which was significantly higher in very high MEPs compared to very low ones regardless the pulse condition (p < 0.001 in SP, p = 0.017 in PP), with the difference being significant also for the high and low trials in SP only.This trend was mirrored in the GMFP for both conditions (p = 0.018, p = 0.008 for SP and PP) and in the LMFP for SP (p = 0.003).In this condition, the TEP mean and area of very high MEPs were significantly stronger than the other labels (p < 0.03).In both pulse conditions, MEP size correlated with the RQS, with the pRQS of higher MEPs being larger than lower ones (p < 0.05, figure 4 fourth row).In SP the pRQS of very high MEPs was also significantly bigger than the pRQS of very low MEPs (p < 0.001).Relationships were also seen in the timefrequency domain: in SP the mean and maximum evoked power was significantly higher in very high compared to very low MEPs trials (p < 0.01).In PP relationships were seen on ITC, with high MEP trials having a stronger maximum and mean ITC than very low trials (p < 0.05).

Discussion
The present work highlighted the impact of the methodological choice of EEG processing and statistical models on the outcomes of brain-state related effects on TMS-evoked responses.It showed how choices in this matter may be responsible for the heterogeneity observed in the literature and how using trials with neural oscillations may reduce the processing effect.With a well-defined processing, we have then presented the effect of alpha oscillations on MEP amplitude and various TEP features in young healthy subjects.Both phase and power modulated MEP amplitude, with phase similarly modulating MEP latency and being robust to experimental changes (i.e.SP vs. PP).The alpha phase also explained some variance of TEP features in the time-domain, whereas power correlated with time-frequency domain features.Again, the pulse condition of TMS (SP or PP) did not result in strong differences on the brain-state-related effects.Finally, strong relationships between MEP size and TEP features of both domains were also observed.

Key role of EEG preprocessing in the assessment of ongoing brain-state
To the authors' knowledge, this is the first time that outcomes based on EEG processing steps have been assessed according to the presence of neural oscillations in specific frequency bands.We observed a higher similarity in results (i.e.phase and PSD estimation) across different pipelines when only trials with neural oscillations were used compared to using all trials or only those without neural oscillations (supplementary figures 2 and 3).Given that the presence of neural oscillations is determined according to the presence of a peak in the PSD spectrum in the frequency of interest when periodic and aperiodic signals are separated [6], the higher correlation between outcomes across processing pipelines may derive from a higher signal-to-noise ratio (SNR) in these trials.This is in line with previous studies showing the importance of PSD and SNR in phase determination [12,56,61], which suggested for instance to measure the phase with multiple filters and looking at the robustness of results [61,62].Separating periodic and aperiodic components of the power spectra is becoming a new trend in the EEG field [17,35,43,63], and our results confirmed its importance in the context of brain-state-dependent stimulation.
However, focusing on trials with neural oscillations does not remove the relevance of specific parameters of EEG preprocessing.When evaluating the effect of brain-state in the alpha band on MEP peakto-peak amplitude in young adults during SP TMS, we obtained different results when comparing the three EEG processing pipelines.Specifically, in addition to our pipeline, we tried to reproduce results of Hussain et al [30] who performed an offline analysis and Madsen et al [41] who had run an online brain-state dependent experiment.The former had observed a significant effect of phase X PSD interaction, with PSD positively correlating with MEP size only in trough trials.The second study only observed an effect of PSD, negatively correlating with MEP amplitude, when only high PSD was considered.With our dataset, we reproduced the results of [30] (figure 3(c)); moreover we observed a main effect of phase, with trough trials eliciting bigger MEPs-this effect was not observed directly by [30], but by other studies performed by the same research group [37,39,40].A main positive effect of PSD was also seen, which was observed by the same research group in another study [38].Results from [41] were not reproduced, but we found a main effect of phase on the MEP amplitude (figure 3(b)).It must be pointed out that the different pipelines used different epoching windows and PSD estimation methods, leading to additional discrepancies in the identification of neural oscillations for the same trial.Although the three models shared the effect of phase, the outcomes were not always consistent.This result is of critical importance as it highlights the significance of EEG processing and makes it very difficult to compare previous results if different methodologies were used to extract the independent variables.Drawing strong conclusions related to electrophysiology becomes even more complicated as different pipelines give different labels in terms of phase for the same trial (i.e.around 40% of trials were labelled differently across approaches) and no real ground truth is available.However, with regard to this important difference in labeling, we observe that most discrepancies are of neighboring labels (e.g. if one approach labels a trial as trough, the other may label it as falling or rising) and only 2% of the trials are estimated as opposite phases.It is possibly due to this that despite the mislabeling, we see a similar phase effect across methodologies.Overall, the present work underlines the relevance of offline analyses to make hypotheses before designing and running online experiments.We believe that such hypotheses should be based on both empirical results and theoretical knowledge.Finally, the statistical model, and most importantly possible variables distributions normalizations used may also have a significant effect on the results reported (supplementary figure 1); this aspect was also recently pointed out by [17].Specifically, we have shown how the proposed generalized linear mixed model with a gamma distribution is better suited for this type of analyses, where symmetric normal distributions violate the physical limits of MEP amplitude that cannot be a negative value.
We have shown that studying the brain-state effect on MEP and TEP features variance is not easily generalizable because of EEG processing and statistics used, but it can also depend on hardware and software differences as well as individual (e.g.anatomical) traits.We encourage using a pragmatic approach for future studies to find a way to maximize the wanted effects.This will be critical for clinical applications and protocols aiming at inducing plasticity.

MEP variance explained by alpha oscillations
Many studies have found an effect of brain state on MEP peak-to-peak amplitude, though with diverging trends (see differences across table 1).Setting aside the effects of processing and confident in the accuracy of the method chosen for the phase estimation, we have focused on the possible explanations of the observed effect.In our analyses we have seen an effect in terms of PSD, which had a positive trend on MEP amplitude in SP; specifically, we observed a clear positive relationship between MEP size and PSD for the trough trials, suggesting that trough trials elicit higher MEPs in presence of high power.The same was seen by [34] and implicitly by [39], which preselected subjects on the basis of power in the alpha band.Focusing on the phase, our results together with the outcomes of figure 3 confirm a literature trend stating that the oscillatory phase of the alpha band in the sensorimotor cortex can explain excitatory states of the corticospinal tract [30,35,39,43].Trough and rising phases are the most excitable and significantly different from the peak and falling phases; this result can be interpreted following the pulsed-inhibition theory, where the trough phase of alpha is excitatory, and the peak phase is inhibitory [40,64].Here, we expand the trough phase to the rising one [35,43].Importantly, the number of trials per phase label is comparable (i.e. between 23% and 27% of trials per label, out of four labels), suggesting that the presence of the different states is similar and that there is not one phase lasting longer than others, which would be independent of the asymmetry of the arc-shaped mu rhythm [65].This result can also be due to the filtering step applied, where the signal is forced to be a pure sinusoid, given the narrow passband bandwidth.The effect of brain-state on MEP amplitude was also evaluated during PP with a SICI paradigm.With SP the corticospinal excitability can be directly assessed; on the other hand, PP stimulation, in addition to triggering the corticospinal tract, pre-activates the intracortical GABAerigc system, which results in a reduced excitability and thus a decrease in MEP amplitude.During PP, we did not reproduce all the effects of SP.This may be due to the decreased intra-and intersubject MEP amplitude variability during PP, as suggested by [40].However, the phase effect remained with peak trials triggering smaller MEPs compared to the falling and rising phase trials, which is a similar effect to that observed in SP.The loss of the PSD main and interaction effects, but not that of phase might give further insights on the physiology of GABAergic inhibition, here probed with SICI.The activation of inhibitory interneurons with the conditioning stimulus reduces the overall corticospinal excitability, but it does not remove the effect of the cortical neurons' excitability state (i.e.phase), which can explain part of the observed variance in MEPs.The GABAergic system seems instead to act upon the PSD effect on corticospinal excitability as during PP the PSD effect is lost.This suggests that the conditioning pulse induces a similar inhibitory effect no matter the size of the neuronal population currently firing.Furthermore, we found a main effect of brainstate phase for MEP latency in SP.Given that the labels showing significant difference are the same seen for the MEP amplitude (i.e.peak vs. trough), it is complicated to state if the shorter latency is correlated to the bigger MEP size, to the alpha oscillation, or a combination of the two.Indeed, MEP size and latency seem to be correlated regardless of the brain state, as also reported in [66].However, in a recent work by Torrecillos et al [32], the authors were able to disentangle the MEP size and the brain-state effect on the latency for the beta oscillations; it remains to be investigated if the physiological model proposed there could be translated to the alpha band.

TEPs features modulation by either alpha phase or power
The literature on the influence of brain-state on TEP features is still very limited compared to the one targeting MEPs (see table 2).Previous studies have often focused on few features, belonging to one general class such as time-domain [15,45] or mean field potential [46]; only Ding et al [13], performed an analysis covering different features domains.However, their protocol was significantly different from ours as they measured brain-state from the occipital cortex, and not from the motor cortex.Overall, we did not reproduce previously seen effects.In our analyses we appreciated the effect of phase in the time domain and of PSD in the time-frequency domain.PSD correlated with evoked power, no matter the frequency band suggesting that if a stronger power is present at baseline, this will be also present after the TMS pulse.Consequently, this will not modulate TEPs dynamics, as confirmed by the lack of effect of PSD in the analysis of RQS.The RQS is a rather novel feature in TEP analyses, and we believe it can be very in studying both the dynamics of the cortical response, with the uRQS, as well as the cortical reactivity with the pRQS.In particular, the higher the pRQS, the more reactive and stable is the cortical response for that specific label (e.g.peak phase), revealing a higher excitability of the targeted neural population [57,58].The higher the uRQS the more similar are the dynamics of the responses evoked with two different labels (e.g. during peak and trough), indicating that the underling physiological characteristics of the stimulated cortical sites are similar (mainly: similar activated neural populations, similar connectivity, similar balance).We also observed higher ITC in high power trials suggesting that when the alpha power is high, the response to TMS can be more consistent across trials.
The alpha phase did not impact the same features, indicating that phase and power have different influences on cortical excitability and dynamics.
Focusing on phase, we saw stronger LMFP and P60 in the trough compared to peak trials (figures 4(a) and (b)), in line with the results [15].This modulation within the early components' relative amplitudes led to significant in evoked dynamics, as revealed by our RQS analysis (figure 4(d)).Considering that (1) P60 has been related to the sensorimotor afferent feedback [16], and that (2) in trough trials we obtained higher MEPs, it is difficult to say which feature is causing the other.Indeed, we have also seen that in presence of higher MEPs, we have higher GMFP, LMFP P60.The same holds true for the time-frequency domain: higher PSD leads to higher MEPs and to higher cortical evoked power, but also higher MEPs correlate with higher evoked power ITC.Importantly, these correlations are present both in SP and PP.Previous comparisons in TEP features between SP and PP, suggest that earlier responses (i.e.20-80 ms) are not affected by the pulse paradigm, whereas later responses (e.g.N100, P180) are significantly reduced in PP [67].
The disentanglement between brain-state and sensorimotor feedback effects on TEP features was beyond the scope of the present study.However, we propose that further analyses on the topic could strengthen the claim with the design of a valid physiological model or a different dataset, using peripheral evoked potential, or with subthreshold intensities.Indeed, we could look for another quality of brain-state that we are not yet controlling, whether a frequency band, a more general microstate, or a different source from where to extract the brain-state (e.g.source localisation) [45].Future studies may try to investigate other frequency bands (e.g.beta oscillations) or to exploit machine learning algorithms to extract potential features.

Conclusion
This work sheds further light on two main aspects related to brain-state-dependent TMS: the effect of EEG processing and the actual effect of brain-state on MEP and TEP features in young healthy adults.Five main points can be drawn: (1) EEG processing is key and minor changes in pipelines steps can lead to major differences in outcomes; (2) choosing only trials with neural oscillations in the frequency band of interest may help diminish variance and as such the impact of different EEG processing methods better allowing to hypothesise neurophysiological mechanisms; (3) alpha phase and PSD as descriptors of brain-state can explain part of the observed variance in MEP and TEP features; (4) experimental manipulation (e.g.PP vs. SP) strongly affects MEP and TEP features at the whole level; however it often does not remove the underlying brain-state-related effect; (5) further investigations are necessary to disentangle the interaction between brain-state, cortical response, corticospinal response, and sensorimotor feedback.Overall, we believe that brain-state dependent stimulation should be further explored, possibly by looking into other brain activity features (e.g.other frequency bands or other regions of interest).For any next step, we suggest benchmarking hypotheses with offline analyses prior to moving to online brain-statedependent experiments-this would be a key step in view of therapeutic experiments with pathological cohorts.
overlapping.Amplitude as square of FFT.Reference for ERD is PSD in 3 s before cue PHASE = NA PSD = negative relationship in SP and SICI (increase of MEP amplitude with ERD).Not in ICF PH x PSD = ANOVA with Bonferroni post-hoc Schaworonkow et al 2018 [37] healthy, right-handed pre-selected in rMT and alpha-peak N = 18 (F = positive negative peak-to-peak amplitude Wilcoxon signed-rank test Thies et al 2018 [38] healthy, right-handed pre-selected in rMT and alpha-peak, N = 16 (F = rmANOVA with post-hoc t-tests correlation per subject tested with one-sample t-test (Continued.) paired t-test in SP rmANOVA with rTMS Bergmann et al 2019 [40] healthy, right-handed pre-selected in rMT and alpha-peak, N = 23 (F = random phase; highest quintile PSD and peak, trough, rising, falling phases) ACTIVITY = rest, eyes open TMS = NR monophasic SP and SICI (ISI = 2 ms) trials, TI = to obtain MEP of 1 mV (50% for condition stimulus), ISI = depends EEG = 64-ch, EMG = rFDI WINDOW: 512 ms SPATIAL F: C3-Hjort FREQ F: two-pass (zero-phase) FIR order 128, 2 Hz around individual alpha frequency PHASE EST: 64 ms edge removal + autoregressive model Yule-Walker order 30 + Hilbert transform PSD EST: Hanning-windowed FFT PHASE = MEP maximal in trough and rising and minimal in peak and falling in SP.No effect in et al 2019 Experiment 3[31] healthy, right-and left-handed, 12 Hz) PHASE EST: 64 ms edge removal + autoregressive model Yule-Walker order 30 + Hilbert transform PSD EST: autoregressive model Yule-Walker order 200 OTHER: demean PHASE = IO-curve shifts to the left for trough compared to peaks; difference between peaks higher at lower stimulation intensities (around 100%rMT) and smaller at higher intensity (around 138%rMT).It peaks at 50% max IO single-pulse, SICI = short-intracortical inhibition; ICF = intracortical facilitation; TI = test intensity, rMT = resting motor threshold, ISI = inter-trial stimulus, reported in seconds; EMG related muscles are reported with r = right, l = left and the abbreviation of the muscle (FDI = first dorsal intraosseus, APB = abductor pollicis brevis).Additional abbreviations: TESA = TMS-EEG signal analyser [44], ERSP = evoked related spectral perturbation, GMFP = global mean field potential, ICA = independent component analysis, ITPC =inter-trial phase coherence.et al 2019[15] healthy, right-handed pre-selected in rMT and alpha-peak, Python 3.8 with scipy, numpy and mne as main packages[54].

Figure 2 .
Figure2.TMS evoked potential features.Top: TMS evoked potentials (TEPs), with the average highlighted in black.Known peaks are pointed at, and their amplitude and latency can be studied.Bottom left: time-domain features are reported such as the global and local mean field potential (GMFP, LMFP).For the latter, in addition to the area under the curve (highlighted in blue), the number of deflections is also evaluated (yellow lines).Bottom middle: regression quality score.The regression quality score look how similar single-trial TEP are to a reference TEP.The score may be a paired RQS if single-trials and reference TEP belong to the same class (top figure) or unpaired RQS if single-trials and reference TEP do not belong to the same class (bottom figure) On the bottom right, time-frequency analyses graphs are reported: at the top evoked power is shown, while at the bottom inter-trial coherence (ITC) is displayed.The full list of measured features is reported in table3.

Figure 3 .
Figure3.Brain-state effect on MEP peak-to-peak amplitude in single-pulse condition according to different approaches.The figure shows the results from the generalized linear mixed models for the alpha band (see tables S1, S5 and S7).On the left the main effect of phase is shown and, on the right, the main effect of PSD (black curve) and its interaction with phase (colored curves).The MEP amplitude used for the models was kept raw, but for consistency and easier view it was z-score normalized per subject for the plot.The PSD used in the model was z-score normalized per subject.(A) Lausanne approach-full statistical results are reported in table S1; (B) Copenhagen approach[41]-full statistical results are reported in table S7; (C) Tübingen approach[30]-full statistical results are reported in table S5.Significant differences are reported.If the main effect of PSD was observed, it is reported, together with its trend at the top left corner of the phase-PSD interaction plot.The estimate of the model is also reported.Color legend: trough phase: green, rising phase: orange, peak phase: violet, falling phase: pink.Significance abbreviations: p < 0.1: '.' , p < 0.05: ' * ' , p < 0.01: ' * * ' , p < 0.001: ' * * * ' .

Figure 4 .
Figure 4. Effects of brain-state and MEP size on TEP features from different classes for different stimulation conditions.Effect of phase is reported on the left column (color legend: trough = green, rising = orange, peak = violet, falling = pink); effect of power is reported in the middle column (color legend: very low PSD = dark blue, low PSD = light blue, high PSD = yellow, very high PSD = orange); effect of MEP size is reported in the right column column (color legend: very low MEP = dark green, low MEP = light green, high MEP = light violet, very high PSD = dark violet).(a) Amplitude of the P60.(b) Local mean field potential (LMFP).(c) Mean evoked power.(d) Regression quality score (RQS); significant differences between paired RQS (pRQS) and unpaired RQS (uRQS) are reported next to the pRQS bar, i.e. the bar obtained when single-trials are regressed over the TEP they belong to-they are represented with the same color of the reference TEP group.The higher the RQS, the more similar evoked dynamics are between the reference TEP and single-trials.Stars indicating significance are colored to identify the uRQS of interest.In these plots, dots represent the mean and the bars the standard errors.Abbreviations: SP = single pulse, PP = paired-pulse, T = trough, R = rising, P = peak, F = falling; vL = very low, L = low, H = high, vH = very high.Significance abbreviations: p < 0.1: '.' , p < 0.05: ' * ' , p < 0.01: ' * * ' , p < 0.001: ' * * * ' from linear mixed effect model.

Table 2 .
Synoptic table of literature review regarding articles where TEP features have been considered in relation to the brain state at the time of TMS stimulation.Both offline and online studies are reported; the type of protocol is highlighted in bold in the column 'Experiment condition' .In this table we report in a summary the following information: reference paper and subjects demographics (i.e.sample size, age, and additional information); experiment condition (i.e.protocol with information of activity performed by subject, TMS related parameters, EEG and EMG electrodes used and targeted cortex area); EEG processing pipelines including both preprocessing, choices of hyperparameters (e.g.window duration of epoch) and how PSD and phase have been estimated) to extract the brain-state; EEG processing (PREP) for TEP features (FEATURES) extraction; results in the alpha band to related TEP features in terms of phase (PH), PSD or interaction effect (PH x PSD); statistics used.If an analysis was not performed by the article, results are reported as NA (not applicable); if a method was not specified, it is reported as NR (not reported).Abbreviations and units used going from left to right and top-bottom: N is the sample size, F is the number of female; age is reported in years, either mean ± SD (range) if available; SP =