Investigating the Origin of TMS-evoked Brain Potentials Using Topographic Analysis

The combination of transcranial magnetic stimulation (TMS) and electroencephalography (EEG) represents an increasingly popular tool to non-invasively probe cortical excitability in humans. TMS-evoked brain potentials (TEPs) are composed of successive components reflecting the propagation of activity from the site of stimulation, thereby providing information on the state of brain networks. However, TMS also generates peripherally evoked sensory activity which contributes to TEP waveforms and hinders their interpretation. In the present study, we examined whether topographic analysis of TEPs elicited by stimulation of two distinct cortical targets can disentangle confounding signals from the genuine TMS-evoked cortical response. In 20 healthy subjects, TEPs were evoked by stimulation of the left primary motor cortex (M1) and the left angular gyrus (AG). Topographic dissimilarity analysis and microstate analysis were used to identify target-specific TEP components. Furthermore, we explored the contribution of cortico-spinal activation by comparing TEPs elicited by stimulation below and above the threshold to evoke motor responses. We observed topographic dissimilarity between M1 and AG TEPs until approximately 80 ms post-stimulus and identified early TEP components that likely reflect specific TMS-evoked activity. Later components peaking at 100 and 180 ms were similar in both datasets and attributed to sensory-evoked activity. Analysis of sub- and supra-threshold M1 TEPs revealed a component at 17 ms that possibly reflects the cortico-spinal output of the stimulated area. Moreover, supra-threshold M1 activation influenced the topography of almost all later components. Together, our results demonstrate the utility of topographic analysis for the evaluation and interpretation of TMS-evoked EEG responses.


Introduction
The combination of transcranial magnetic stimulation (TMS) and electroencephalography (EEG) is becoming an established tool in experimental and clinical brain research, Communicated by Carlo Miniussi. Dominika Sulcova dominika.sulcova@uclouvain.be 1 provide information about the strength of the TMS-evoked cortical response, the location of the sources generating this response and, consequently, the functional connectivity across the TMS-activated cortical network (Bortoletto et al. 2015;Rogasch and Fitzgerald 2013).
However, several issues must be tackled in order to decipher the information held by TEPs. First of all, it is now widely recognized that TEPs reflect several parallel processes of heterogenous origin. In addition to the genuine TMS-evoked cortical activity, they also include somatosensory and/or auditory-evoked activity related to the fact that TMS pulses unavoidably activate a range of sensory receptors (Biabani et al. 2018;Conde et al. 2019;Rocchi et al. 2020). This presumed peripheral contribution to TEPs received much attention in recent years, culminating with a paper of Conde et al. (Conde et al. 2019) that compared brain responses elicited by active TMS and realistic sham stimulation. The group reported high similarity of both evoked potentials, suggesting that TEPs predominantly reflect sensory-evoked brain activity. Although their findings have been widely disputed (Belardinelli et al. 2019) and multiple research groups have since striven to prove that TEPs do contain measurable TMS-evoked cortical activity (Ahn and Fröhlich 2021;Freedberg et al. 2020;Rocchi et al. 2021), it still indicates that the functional significance of each identified TEP component must be carefully considered.
The situation further complexifies when it comes to TEPs evoked by stimulation of the primary motor cortex (M1). In any other cortical area, TMS leads to the activation of corresponding cortico-cortical and cortico-thalamic circuits. While this is also true for sub-threshold M1 stimulation, supra-threshold stimulation of the area causes an additional excitation of the cortico-spinal tract that provokes a contraction in the target muscle. This peripheral event can be easily registered as a motor evoked potential (MEP), which represents a readily available measure of the overall motor pathway excitability. Unsurprisingly, this outstanding property of M1 has made it the most popular cortical target for TMS. However, the muscular contraction evoked by M1 stimulation elicits an additional time-locked, re-afferent somatosensory brain response that is mixed in the resulting TEP (Gordon et al. 2018;Komssi et al. 2004), thus complicating its interpretation.
Finally, assuming that one successfully identifies a set of TEP components that can be linked to the cortical process of interest, yet another challenge arises during analysis. Like any other type of event-related brain potentials (ERPs), TEPs are most often evaluated in the time domain where the amplitude and latency of individual peaks represent the principal measured variables. Summarizing the TEP waveform as a series of peaks offers the advantage of substantially reducing the number of considered variables, but it usually requires prior knowledge on the time course and topographic distribution of the elicited responses to allow a preliminary selection of channels and time windows of interest.
Another way to evaluate ERPs, proposed by Murray et al. (2008), is to consider the temporal evolution of scalp topography of the evoked response rather than its amplitude. This approach is based on the notion that the potentials recorded at the scalp reflect the momentary state of global neuronal activity, with all simultaneously active sources contributing with different weights to the signal recorded by each electrode (Vaughan 1982). It implies that two different topographies recorded from the same subject represent the activation of two distinct sets of cortical sources (Michel et al. 2004). Therefore, by searching for topographic differences across time and experimental conditions, it is possible to identify time intervals of differential cortical activity in a purely data-driven fashion. Topography can be also compared across timepoints within a single dataset in order to identify clusters corresponding to the period of activity of a certain set of sources. This approach, known as microstate analysis (for review see (Khanna et al. 2015;Michel and Koenig 2018), stems from an early study (Lehmann et al. 1987) that identified short intervals of quasi-stable topographies -microstates -in the resting-state EEG. It was later suggested that microstates might represent the building blocks of information processing (Lehmann 1990). In evoked potentials, microstate segments are closely aligned to observed ERP components and are thought to correspond to a sequence of discrete processing steps triggered by the stimulation (Murray et al. 2008).
To our knowledge, neither analysis of topographic similarity nor microstate analysis are commonly used to characterize TEPs, despite several non-negligible strengths of these methods. By including the whole set of channels in the form of a single multidimensional vector, a statisticsbased topographic analysis avoids the necessity to arbitrarily choose regions of interest. Moreover, by considering the voltage landscape as a whole, this approach effectively deals with the bias introduced by the choice of the reference electrode. The placement of reference inevitably influences the scalp voltage amplitude (Geselowitz 1998;Lehmann and Skrandies 1980) and by consequence the character of differences observed between tested conditions, as elegantly illustrated by Murray et al. (Murray et al. 2008).
TMS pulses delivered over two different cortical targets can be expected to elicit cortical responses that differ at least partly in terms of activated cortical networks. However, the stimulation is also likely to evoke a set of unspecific cortical responses due to the activation of highly distributed connections, as well as the concomitant triggering of somatosensory-and auditory-evoked responses. Furthermore, varying the intensity of the TMS pulse can be anticipated to differentially modulate the magnitude of these different response components, and this effect might be especially remarkable when stimulating M1, depending on whether the stimulus activates the cortico-spinal tract. In the present study, we tested whether the topographic analysis of TEPs elicited by stimulation of two distinct cortical targets, namely M1 and the angular gyrus (AG), or using two distinct stimulation intensities (M1 stimulation below and above the threshold to elicit MEPs) can be used to disentangle target-specific and target-unspecific components of TEPs and, thereby, gain information on the properties of the cortical networks activated by TMS.

Ethical approval
All experiments were conducted according to the latest revision of the Declaration of Helsinki. The protocol and all procedures were approved by the local Ethical Committee (Commission d'Ethique Biomedicale Hospitalo-Facultaire). All participants were thoroughly informed of the protocol and applied procedure, gave a written informed consent before taking part in any experiment, and were financially compensated for their participation in the study.
Participants 20 healthy right-handed volunteers (10 males, mean age 24) took part in the study. All participants filled a detailed questionnaire and undertook a brief neurological examination to exclude candidates with any contraindication for experimental procedures, such as epilepsy (or familiar history thereof), presence of metal in the body (e.g. insulin pump, pacemaker, metal prothesis), known neurological pathology relevant to the study and use of medication alternating cortical excitability (Rossi et al. 2009).

Experimental design
The data analysed in the present study corresponds to two baseline datasets acquired in a larger study assessing the pharmacological modulation of TMS-evoked potentials (TEPs). The baseline datasets were acquired in two sessions separated by at least one week, conducted in the morning while keeping onset times as similar as possible for each subject. The subjects were asked to sleep sufficiently and refrain from alcohol consumption on the night prior to the experiment, and to avoid caffeinated beverages in the morning before the experiment.
TEPs were recorded following TMS stimulation of the left (dominant) M1 and the left AG. TMS pulses delivered over M1 were delivered using both a sub-threshold and a supra-threshold stimulation intensity, the latter evoking clear motor contraction in the target muscle, monitored using a surface EMG recording.

EEG recording
The EEG was recorded using a NeurOne EEG system (Bittium NeurOne Tesla; Bittium Corporation, Oulu, Finland) and a 32 channel EEG cap mounted with TMS-compatible Ag/AgCl electrodes (EasyCap GmbH, Herrsching, Germany). The signal was referenced to both mastoids and recorded at a 20 kHz sampling rate with a 5 kHz low-pass filter and DC filter. Electrode impedances were kept below 5kΩ.
A thin layer of plastic was placed over the electrodes and subjects wore an extra rubber cap on top to minimize artifacts caused by the direct contact of electrodes with the TMS coil, and to reduce the bone conduction of the vibration generated by the TMS stimulation. In addition, during the EEG recordings, subjects listened to a masking noise created from an original recording of the TMS click to further reduce possible auditory artifacts.

Neuronavigation and TMS
A 3D T1-weighted structural magnetic resonance image (MRI) of each participant's whole brain was acquired beforehand at the Department of radiology of the Cliniques universitaires Saint-Luc (1 × 1 × 1 mm; 3 T Achieva; Philips Healthcare, Amsterdam, The Netherlands). A Visor2 neuronavigation system (Visor 2.3.3; Advanced Neuro Technologies, Enschede, The Netherlands) was used to verify the accurate placement of the TMS coil and to monitor its position throughout the experiment. A 3D model of the scalp and brain surface was reconstructed based on the individual MRI data and co-registered with the real space using landmark-based markers (nasion and tragi) and head-shape matching (Gugino et al. 2001). The positions of the subject's head and TMS coil were continually registered with a Polaris infrared optical tracking system (Polaris Spectra; Northern Digital Inc. Europe, Radolfzell, Germany).
Both cortical targets (M1 and AG) were identified on subject's brain model and their location was saved in the neuronavigation system to keep the stimulation sites identical across experimental sessions. The TMS pulses were delivered using a figure-of-eight TMS coil having an outer diameter of 75 mm (C-B60; MagVenture, Farum, Denmark), tangentially placed on the scalp and centred over the stimulation target. The stimulation consisted of biphasic complete removal of the artifact in the signal at most channels, however at some we could still observe a leftover tail of the TMS-evoked muscle artifact with a sharp front edge that could pose a problem for the frequency filter. For that reason, an Independent Component Analysis (ICA, FastICA algorithm (Hyvärinen and Oja 2000)) was performed to remove components containing this sharp tail, and in this step epochs of all categories from one recording were concatenated and treated together. Data were then bandpass filtered between 0.1 and 80 Hz (Butterworth, 4th order) and notch filtered at 50 Hz (FFT linear filter, notch width 2 Hz, slope cut-off width 2 Hz). In the next step, individual epochs were visually inspected to remove those containing excessive noise or isolated outbursts thereof. Any trial associated with a missed TMS stimulation (in case of an unexpected head movement causing the coil to slip etc.) was also discarded. On average, the final number of remaining epochs per stimulus category was 48 ± 7 (SD) and there were no significant differences between stimulation conditions.
A second round of ICA was then used to remove remaining non-neural artifacts such as eye blinks, horizontal eye movements, tonic muscle activity and electrode noise. Artefactual components were first identified automatically with the Multiple Artifact Rejection Algorithm (MARA), a linear classifier trained on adult EEG data (Winkler et al. 2011). The pre-selected components were evaluated based on their topography, time course, and power spectrum density (PSD) to remove only components with clearly artefactual origin. Eye blinks and horizontal eye movements were identified by their specific frontal topography, irregular onset and steep decay of the PSD with a high content of slow oscillations. Tonic muscular contraction was easily distinguished by its marginal topography, no remarkable time-locked activity and a high content of high frequency oscillations. Electrode noise was identified based on topography limited to a single channel, often associated with sudden high-amplitude voltage changes. This second ICA was also used to remove leftovers of the TMS artifact caused by TMS-induced contraction of muscles under the stimulation coil. These artifacts usually showed a bipolar topography in the proximity of ipsilateral temporal or neck muscles.
Since M1 and AG TEPs were recorded within the same stimulation blocks, the majority of above mentioned extra-cerebral artifacts were shared across both datasets. However, certain artifacts, such as leftover activity due to muscle contraction, were often specific for the stimulation site. For that reason, an ICA matrix was calculated for each target and the datasets were treated separately. To minimize the risk of introducing differences due to variations in the removal of these artifacts, we verified that the same number and type of common artifacts were removed from M1 and AG TEPs. It is also necessary to emphasize that we took pulses delivered manually using a MagPro X100 magnetic stimulator (MagPro X100 with MagOption; MagVenture, Farum, Denmark), inducing electrical current with anteriorposterior posterior-anterior (AP-PA) direction.
M1 stimulation. The stimulation target over the left M1 was defined as the most optimal position to elicit motorevoked potentials (MEPs) in the right first dorsal interosseous muscle (FDI). The optimal coil orientation was adjusted individually, with the handle always pointing laterally and posteriorly with an approximate 45° angle from the midline. The resting motor threshold (rMT) was set to the minimal stimulation intensity eliciting MEPs with an amplitude of least 50 µV in at least 5 trials out of 10 (Conforto et al. 2004). The mean rMT (in % of maximum stimulator output ± SEM) was 49.3 ± 1.4% and 49.8 ± 1.4% in the first and second experimental sessions, respectively. TMS over M1 was delivered in three blocks of 60 stimuli with a random 4-6 s inter-trial interval while simultaneously recording EEG and EMG. Three categories of stimuli were applied in a semi-random order, with a total of 60 stimuli of each category (20 per block): (1) supra-threshold single TMS pulses at the intensity of 120% rMT, (2) sub-threshold single TMS pulse at the intensity of 80% rMT and (3) paired TMS pulses combining the two types of TMS stimuli. In this study, only TEPs elicited by single supra-threshold and sub-threshold TMS pulses were analysed.
AG stimulation. The stimulation target over the left AG was identified in the individual MRI based on the position of anatomical landmarks, including the Superior temporal sulcus and the Intraparietal sulcus. The coil was placed with the handle pointing downwards and approximately 10° towards the midline, in a position that corresponds to the orientation of the coil axis across the superior temporal sulcus and perpendicular to surrounding gyri. TMS was applied in a single block of 60 stimuli delivered at 100% rMT with a variable inter-trial interval (4-6s).

Pre-processing of the EEG data
The continuous EEG recordings were processed offline using Matlab (MathWorks, Inc., Natick, Massachusetts, United States), Letswave 6 (an open-source EEG/EMG signal processing toolbox, https://www.letswave.org/) and custom scripts. We generally followed the pipeline introduced by Rogasch (Rogasch et al. 2014). The EEG signals were first re-referenced to the average of all scalp channels, the DC shift was removed and a linear detrend was applied. The data were cut into epochs from − 1000 ms to 2000 ms relative to each TMS pulse. The large artifact caused by the TMS was removed and replaced using cubic interpolation between − 5 and + 10 ms, and the signal was down-sampled to the sampling rate of 2 kHz. The interpolation led to a

Topographic analysis of TEPs
All steps of the topographic analysis were performed in Ragu, a free Matlab-based toolbox for the analysis of ERPs (Koenig et al. 2011). Ragu uses a multi-variate approach that considers the entire scalp voltage distribution and allows to make inferences about group differences based on randomization statistics. Our analysis pipeline followed the steps indicated in (Habermann et al. 2018).
Test for outliers. Prior to any further analysis, compared datasets were tested for outliers. All subjects were represented in the form of a single vector composed of all available datapoints. Subjects with datapoints unlikely falling in a normal distribution were automatically identified based on the Mahalanobis distance between subjects (Wilks 1963) and excluded from further analysis.
Test for topographic consistency. A Topographic Consistency Test (TCT) of the TEPs obtained for each session and each type of stimulus, and the obtained results were then combined across sessions. The TCT identifies time intervals that show strong topographic homogeneity across subjects. Changes in the TEP waveforms that are observed in these intervals can thus be reliably linked to an effect of the experimental manipulation on the TMS-evoked activity. Conversely, the TCT helps to avoid attributing any meaning to changes observed at times of excessive inter-individual topographic variability. The TCT implemented in Ragu uses the GFP of the group-level average TEP data as the quantifier, which is then tested for significance against an empirically obtained distribution of GFP values corresponding to the null hypothesis. This distribution was computed at each timepoint using 5000 permutations, with each iteration randomly scrambling the individual-level signal across all electrodes to create random topographies with a constant GFP and recalculating the group average.
Analysis of topographic similarity. To assess the degree of similarity between topographies of two TMS stimulation conditions (AG vs. sub-threshold M1 stimulation; subthreshold vs. supra-threshold M1 stimulation), we used the measure of Global Dissimilarity (DISS, Eq. 2) also introduced by (Lehmann and Skrandies 1980). The DISS corresponds to the distance between two multidimensional vectors each representing one entire scalp voltage distribution across electrodes. The DISS value ranges between 0 (perfect topographic homogeneity) and 2 (complete topographic inversion). For each of the planned comparisons, a single value of DISS was obtained for each timepoint from the group-level averaged TEPs of the compared conditions. A point-by-point statistical randomization test was conducted to determine the time windows of significant difference between conditions. In Ragu, this procedure is referred to as a Topographic Analysis of Variance (TANOVA), even great care not to discard components that were suspected to possibly include TMS-evoked cortical activity, such as topographically widespread changes in signal in the background of a localized artifact, or any component that displayed non-explainable time-locked activity in the average time course. Thus, the cleaned EEG signals most probably included some weak remaining artefactual signals.
As a last step, the EEG epochs were baseline corrected to the interval from 200 ms to 5 ms before the TMS pulse and averaged for each subject and condition.

TEP peak labelling
The TEP waveforms were evaluated within the time interval from 10 to 300 ms post-stimulus. For each stimulation type, individual components were first identified as the local maxima in the time course of the Global Field Power (GFP) that was calculated from the grand average TEPs. GFP corresponds to the standard deviation of voltage across all channels that can be plotted as a function of time (Eq. 1). As such, it conveys information about the overall strength of time-locked TMS-evoked responses, which itself corresponds to the degree of synchronisation of TMS-evoked cortical activity (Lehmann and Skrandies 1980). The latencies of GFP peaks were extracted automatically and used to search for corresponding deflections in the signal recorded at the vertex electrode. TEP components were then named according to the polarity and approximate latency of these deflections, while keeping in line with the nomenclature established in the literature when possible.

Planned comparisons
To evaluate the effect of stimulation site, TEPs elicited by AG stimulation were compared to the TEPs elicited by subthreshold M1 stimulation. This paired comparison was justified by the fact that sub-threshold TMS over M1 does not evoke any contraction in the target muscle, and therefore does not trigger EEG responses related to the processing of recurrent somatosensory signals originating from the TMS-evoked contraction of the hand. TEPs elicited by subthreshold M1 stimulation should thus be more comparable to TEPs elicited by stimulation of non-motor cortical areas such as AG.
Then, we compared the TEPs elicited by sub-threshold M1 stimulation to the TEPs elicited by supra-threshold M1 stimulation, with the intention to disentangle the contribution of somatosensory processing to the M1 response.
For both comparisons, the datasets from the two experimental sessions were included and the factor session was added to the statistical analysis to control for possible testretest variability.
per condition for each segment of statistically significant topographic dissimilarity. The voltage difference was compared at each channel and obtained t values were plotted as t-maps. These t-maps allow to visualize where and how much mean topographies differ by comparing the average difference to the variance across observations. Reported p values were obtained by averaging p values provided by TANOVA across all timepoints in each time interval.
Microstate analysis. A microstate analysis of the TEP waveforms was used to assess the temporal properties of observed topographic changes. The optimal number n of microstate classes was determined using the cross-validation method proposed by (Koenig et al. 2014), as illustrated in the lower panel of Fig. 1. First, the average TEP waveforms of individual subjects were randomly split into a training and a testing dataset. Then, the average of training TEP waveforms was used to compute microstate models with the number of classes ranging from 3 to 12, and the mean spatial correlation C i of each model with the average of the testing dataset was computed according to Eq. 3. This process was repeated 50 times for variable subsets of the data. The C i from all runs was averaged and compared across models to indicate the simplest one yielding the C value in the plateau range (the 'elbow position') and the corresponding n of microstate classes was retained for the subsequent microstate labelling and analysis.
The n class maps were identified in a dataset created by concatenating the group-level average TEP waveforms of compared stimulation conditions (this process is illustrated in the upper panel of Fig. 1). Each timepoint in each group-level average TEP waveform was then labelled with the class map that was most correlated with the momentary topography, thus segmenting the TEPs into non-overlapping intervals or microstates. For the purpose of this segmentation, we chose to use the k-means clustering algorithm (Murray et al. 2008) with 250 iterations. In each cycle, the algorithm drew n random template maps (topographies at n different timepoints) from the concatenated dataset and computed the spatial correlation C between each template and all the remaining topographies in the dataset. In the next step, each timepoint in the dataset was labelled with the template map that yielded the highest C, and the Global Explained Variance GEV i (Eq. 4) was computed for the given set of templates. All template maps were subsequently re-defined by averaging topographies over all the timepoints labelled with the same template. C and GEV i were re-calculated and the whole process was repeated until no additional increase in GEV i was observed. Obtained GEV i was then associated with the final, refined set of microstate class maps, and the next iteration of the clustering process began with yet another random sample of n topographies. Finally, the set of n class maps that yielded the highest overall GEV i though it is a non-parametric test unrelated to an analysis of variance (ANOVA). Indeed, the TANOVA computes for each time point an empirical distribution of DISS values under the null hypothesis by running 5000 cycles of randomly assigning individual TEP waveforms to the two compared conditions and re-calculating the DISS. The actual DISS values are then compared at each timepoint to the values from this empirical distribution, thus yielding a time-varying p value.
Finally, the p values obtained for each time point were corrected for multiple comparisons by applying a procedure named "Global Count Statistics", which compares the obtained number of sub-threshold p values with the count of sub-threshold p values in an empirically constructed distribution of p values under the null hypothesis. The level of significance was set to 0.05.
To further characterize observed topographic differences across conditions, mean scalp topographies were calculated Fig. 1 Microstate analysis. The upper panel: microstate labelling using k-means clustering algorithm. This process is applied at the level of grand average datasets. In a single iteration (left), a set of microstate class maps is obtained by randomly selecting template topographies of n timepoints (1), labelling all timepoints with a single template topography that yields the highest spatial correlation (2) and refining those templates by averaging all timepoints with the same label (3). Each refining cycle, the Global Explained Variance of the template set is calculated (GEV i ) and the set is considered final when no increase in GEV i can be observed. Therefore, the outcome of one iteration is a set of refined template maps and an associated value of GEV i . GEV i obtained in all iterations is then used to identify the most fitting set of template maps, that model is retained, and statistics are computed based on the corresponding microstate labelling (right). Lower panel: choosing the optimal number of microstate classes using cross-validation. In a single iteration (left), individual averaged data are split to a training and testing dataset (1), the training dataset is used to prepare a range of microstate models with increasing n (using the process shown in the upper panel) (2) and global spatial correlation C i is computed between each model and the testing dataset. C i obtained in all iterations is averaged and the n of the simplest model yielding C in the plateau range is kept for the final analysis (right) Equations of (1) Global Field Power GFP, (2) Global Dissimilarity DISS, (3) Spatial Correlation C, and (4) Global Explained Variance GEV. GFP, DISS and C are calculated for each timepoint t , whereas GEV averages across all timepoints. n is the number of electrodes, u j and v j represent the voltage measured at one electrode in two different conditions, while − u is the mean voltage computed from all electrodes in a single condition. M t stands for the class map attributed to the timepoint t .

TEPs elicited by supra-threshold M1, sub-threshold M1 and AG stimulation
For each stimulation condition, clear peaks were observed in the TEP waveforms. These peaks differed in their topography and/or latency across stimulation conditions and showed variable topographic consistency.
In the TEP elicited by supra-threshold M1 stimulation, we identified a sequence of five highly consistent peaks that were labelled N17 (peak latency 17 ms), P30 (32 ms), P60 (57 ms), N100 (111 ms), and P180 (194 ms). Associated peak topographies are displayed in Fig. 2a. An additional was retained and backfitted to the group-level average TEP waveforms.
The final microstate model labels every timepoint of the group-level average TEP waveforms with a single microstate class. Using this information, it is possible to compare temporal properties of those classes that are shared across compared datasets. To this end, we considered the following microstate features: (1) mean duration, (2) onset and offset and (3) centre of gravity. To compare these measures between conditions, we quantified the observed effect by computing the difference, in other words the variance, across conditions. Individual datasets were randomized with 5000 permutations to create the distribution of the given measure under the null hypothesis, and the original quantifier was compared to this distribution. The probability of the quantifier being compatible with the null hypothesis was computed as the proportion of trials in the whole empirical distribution that gained values larger than the quantifier.

Analysis of the topographic dissimilarity between sub-threshold M1 TEPs and AG TEPs
Prior the analysis, two outliers were identified and excluded. The DISS analysis identified three intervals of significant topographic dissimilarity between the TEP waveforms elicited by sub-threshold M1 stimulation vs. AG stimulation. Since these time windows partially overlapped with intervals of topographic inconsistency identified in both datasets, the inconsistent intervals of both time windows were excluded from further analyses (Fig. 3a).
The first interval of significant dissimilarity (21.5-26 ms, p < 0.001) partially covered the P30 component of M1 TEPs and the P25 component of AG TEPs. The second interval (54.5-80.5 ms, p < 0.001) included the P60 component of M1 TEPs and the N75 component of AG TEPs. For both intervals, the M1 TEP topography featured a positive peak was identified around 45 ms and labelled N45. The topography of N45 was masked by an overlapping P30-P60 complex but was nevertheless distinguishable as a negative deflection superimposed on the positive peak.
Five peaks were also identified in the TEP waveforms elicited by sub-threshold M1 stimulation: P30 (27 ms), N45 (46 ms), P60 (61 ms), N100 (118 ms) and P180 (184 ms). The TCT showed low topographic consistency between 10 and 21 ms, i.e. at latencies corresponding to the N17 observed following supra-threshold M1 stimulation. Low topographic consistency was also observed 33-50.5 ms, which included the N45. The latencies of these peaks were very similar to the latencies of the peaks elicited by suprathreshold M1 stimulation. However, the peak voltage distribution of all peaks except the P30 visibly differed from that of the corresponding peaks for supra-threshold M1 stimulation, as shown in Fig. 2b.
In the TEP waveforms elicited by AG stimulation, 6 components were identified: P25 (21 ms), N40 (37 ms), P55 (55 ms), N75 (75 ms), N100 (114 ms) and P180 (192 ms). Peak topographies are shown in Fig. 2c. The TCT revealed associated with the class 3 microstate that was also ascribed to the N75 component of AG TEPs described above. The P180 was labelled with a microstate class composed of a positive peak centred over the vertex (Fig. 3b: class 5). The periods of lower GFP between components were segmented into short intervals of transitory microstates.

Analysis of the topographic dissimilarity between sub-threshold and supra-threshold M1 TEPs
One outlier was identified and excluded from further analysis. The topographic dissimilarity analysis revealed that the topographic distribution of sub-threshold vs. supra-threshold M1 TEPs were dissimilar over most of the explored time window, while no differences were detected between TEPs from the two sessions (two short intervals of significance, 77-87 and 130.5-138 ms, did not survive the correction for cluster duration; Fig. 4a). Further exploration of the topographic differences related to stimulation intensity was restricted to the intervals showing topographic consistency across both stimulation conditions. The first identified interval extended 21.5-33 ms. The second period of topographic dissimilarity (51-221.5 ms) encompassed several TEP peaks having very different topographies. Therefore, it was segmented into three separate intervals of interest (interval 2: 51-67, interval 3: 81-153.5 and interval 4: 167-217 ms) by lowering the significance threshold and only including timepoints with a p value lower than 0.005. A fifth interval of significant dissimilarity was identified towards the end of the explored time window (230.5-267.5 ms; Fig. 4a).
The first interval of significance (21.5-33 ms, p < 0.01) covered the P30 component of M1 TEPs. Despite the significant dissimilarity value, the mean maps showed relatively similar topographies with a maximum close to the vertex for both TEPs, shifted towards left central electrodes for the sub-threshold TEP, and more symmetric and widespread for the supra-threshold TEP. This was confirmed by the t-map with negative values over left central electrodes and widespread positive values over the right hemisphere.
The second interval (51-67 ms, p < 0.001) contained the component P60. For sub-threshold stimulation, the mean topography was very similar to that of the first interval. The third interval of dissimilarity (88.5-106.5 ms, p < 0.05) covered the N100 component of both M1 TEPs and AG TEPs. In this interval, the mean topographic maps of M1 and AG TEPs were visually rather similar, both datasets showing maximum positive voltages at posterior electrodes and maximum negative voltages at anterior electrodes, but with a lateralized distribution towards the right for M1 TEPs, and a more symmetrical distribution across the two hemispheres for AG TEPs. The t-map shows a diffuse pattern with the positive t-values over the right hemisphere, and negative t-values over the left hemisphere.
No significant topographic differences were observed between TEPs recorded in the two sessions. The p value obtained by TANOVA for the factor session (Fig. 3a) exceeded 0.05 in two short intervals (149.5-152 and 196.5-197.5 ms) that were not retained when corrected for cluster duration.

Microstate analysis of sub-threshold M1 TEPs and AG TEPs
The microstate analysis yielded six microstate class maps explaining 90.41% of the total variance (Fig. 3b). Up until approximately 80 ms post-stimulus, the M1 TEP and AG TEP waveforms were associated with different microstates.
For the M1 TEP waveform, the P30 and P60 were labelled with the same microstate class having a positive maximum over the left central electrode C1 (Fig. 3b: class  4). The intermediate component N45 was not evaluated due to its low topographic consistency.
For the AG TEP waveform, the P25 was labelled with a microstate class having a bipolar topography negative over left frontal electrodes and positive over posterior electrodes ( Fig. 3b: class 1). The later N75 was characterized by a microstate class showing a widespread frontal negativity maximal over the right hemisphere ( Fig. 3b: class 3). The N40 and P55 were not evaluated due to their low topographic consistency.
The later TEP components, i.e. the N100 and P180 elicited by both M1 and AG stimulation, were composed of the same sequence of microstates with similar latency and duration. Longer intervals of high GFP were associated with one dominant microstate class per component. The N100 was

Microstate analysis of sub-threshold and suprathreshold M1 TEPs
The microstate analysis yielded 8 microstate class maps explaining 91.12% of the variance (Fig. 4b), and confirmed the topographic dissimilarity analysis by attributing different microstates in the compared datasets to the corresponding intervals of high GFP that define TEP components. A single exception was the component P30 that featured the same class map in TEPs evoked by sub-and supra-threshold stimulation. While the supra-threshold TEP waveform was composed of a sequence containing all eight microstate classes (Fig. 4b), the topographically consistent sections of the sub-threshold TEP featured only three prominent microstate classes (Fig. 4b: classes 2, 4, and 7).
In the next step, temporal features of the microstate classes that figure in both datasets were compared between supra-and sub-threshold TEPs. The microstate class 2, which labelled the P30, showed borderline significant earlier onset (p = 0.05) and significantly earlier centre of gravity supra-threshold TEP, it corresponded to a widespread frontal negativity peaking over the left hemisphere. The t-map was characterized by maximally negative t-values over the left central electrode C1 (t = -27.05).
The fourth (167-217 ms, p < 0.001) and fifth (230.5-267.5 ms, p < 0.05) intervals spanned the P180 component. In the earlier window, the topography of the sub-threshold TEP was centred over the vertex whereas the topography of the supra-threshold TEP was widespread slightly shifted towards the right hemisphere. The t-map shows a clearly defined pattern of difference with maximum negative t-values at CP1 (t = -23.12) and the maximum positive t-values at F4 (t = 17.91). In the later window of significance, the mean topographies of both sub-threshold and supra-threshold TEPs were relatively symmetric and maximal over the vertex, but the t-map clearly showed that positive scalp voltages spread more posteriorly in the supra-threshold TEP as compared to the sub-threshold TEP. AG TEPs, the t-map closely mirrored the voltage distribution of mean topographic maps, suggesting a clear separation of M1-evoked vs. AG-evoked cortical responses. By contrast, in the following interval lasting until approximately 80 ms and covering the P60 component of M1 TEPs and the N75 component of AG TEPs, the t-map highly resembled the one obtained in the previous interval but began to differ from current mean topographic maps. This indicates that the set of cortical sources active during the first interval might remain active during this second interval but the signal becomes mixed with additional background activity that might be shared across both datasets. The subsequent microstate analysis labelled the M1 TEP components P30 and P60 and the AG TEP components P25 and N75 as distinct classes, further confirming the spatiotemporal specificity of these early-latency TEP components.
Our results and their interpretation are in line with previous studies using sham TMS stimulation. Herring et al. compared TEPs evoked by stimulation of the visual cortex to the EEG response elicited by sham (active TMS delivered at the shoulder). They found that the evoked potentials became highly similar around 80 ms after stimulation onset (Herring et al. 2015). Using a similar design while stimulating M1, Biabani et al. reported high temporal and spatial correlation between brain responses to M1 stimulation and sham stimulation starting approximately 60 ms post stimulus . More recently, Freedberg et al. compared TEPs across several cortical targets (M1, dorsolateral prefrontal cortex and posterior parietal cortex) and between active and sham stimulation using cosine similarity, which is another amplitude-free measure that allows to evaluate the global shape of stimulus-evoked brain responses (Yaffe et al. 2014). Like us, they found low spatial similarity across TEPs until approximately 60 ms, after which topographic features started to converge, becoming highly similar around 100 ms (Freedberg et al. 2020).
However, it must be pointed out that both datasets showed early periods of topographic inconsistency that were not included in the comparison, thereby limiting our conclusions to a subset of TEP components. We cannot exclude that excessive inter-subject variability in TEP topography indicates an absence of meaningful time-locked activity at these latencies. Or it could correspond to intervals of weak amplitude that is lost in the individual noise, especially when signal averaging relies on a relatively low number of recorded trials (Kerwin et al. 2018), as was the case in the present study.
Finally, it is necessary to acknowledge that early TEP components are highly susceptible to the contamination by artifacts related to the TMS-induced contraction of scalp, face or neck muscles (Rogasch et al. 2014). Targeting different cortical areas requires distinct coil placement and likely (COG, p < 0.05) in the sub-threshold TEP as compared to the supra-threshold TEP. The microstate class 4, which covered the entire duration of N100 in the sub-threshold TEP and only its beginning in the supra-threshold TEP, showed a significantly longer duration (p < 0.05) and later COG (p < 0.05) in the sub-threshold TEP. Lastly, the microstate class 7, which encompassed the entire P180 duration in the sub-threshold TEP, appeared only at its end and was shorter lasting (p < 0.05) in the supra-threshold TEP. For classes 4 and 7, neither onset nor offset latencies were significantly different between sub-and supra-threshold TEPs.

Target specificity of TEPs and the contribution of sensory-evoked responses
In the first part of this study, topographic analysis was used to compare TEPs elicited by stimulation of the primary motor cortex and the angular gyrus, with the aim of distinguishing components corresponding to the genuine cortical activity triggered by TMS stimulation from unspecific components predominantly reflecting cortical activity generated by the peripheral activation of sensory receptors. We assumed that stimulation of each of the two cortical targets would evoke a specific pattern of activation reflecting the direct cortical response to TMS, entangled with additional sensory-evoked activity that might vary in its strength but would have a very similar topographic pattern in both datasets.
Early-latency TEP components. Taken together, our topographic analyses indicate that early-latency TEP components predominantly reflect target-specific activity generated in the proximity of the stimulation site. As such, they likely correspond to the genuine cortical-evoked brain activation that may contain valuable information about the properties of the cortex targeted by TMS and its associated functional networks. While the very early responses, such as the P30 component elicited by M1 stimulation and the P25 component elicited by AG stimulation, appeared to be highly specific for the stimulated region, our data also suggests that other sources of cortical activity, possibly corresponding to early-latency components of a somatosensory-evoked response, might be entangled in the signal of interest from approximately 50 ms. Therefore, any later TEP components should be interpreted with caution. Yet, our conclusion clearly argues against the notion that all TEP components are to a great extent generated by unspecific sensory-evoked activity (Conde et al. 2019).
In our dataset, during the earliest time interval of topographic dissimilarity identified by TANOVA encompassing the P30 component of M1 TEPs and the P25 component of reported that efficient masking allowed to completely suppress the contribution of sensory-evoked potentials to the recorded TEPs, a statement that is problematic for several reasons. First, because the efficiency of auditory masking varies across subjects and the coverage might be insufficient in individuals requiring higher intensities of TMS because of a high threshold to elicit motor-evoked potentials (Ter Braack et al. 2015). Second, there is currently no efficient mean to suppress sensations produced by direct stimulation of peripheral afferents innervating the scalp. In a recent study, Ross et al. attempted to bypass this issue by applying an independent component analysis (ICA) to TEPs evoked by unmasked supra-threshold stimulation of M1 and AG in order to isolate and remove signals attributed to auditoryevoked responses (Ross et al. 2022). Using this approach, they recovered TMS-evoked cortical activity embedded in the auditory-evoked potential, revealing a new component peaking approximately 120 ms post stimulus and showing similar topography across stimulated cortical regions. Nevertheless, it remains uncertain to which extent ICA is able to separate cortical-and sensory-evoked activity, because the assumption of temporal independence of contributing sources is likely violated due to the time-locked character of all interlaced signals ).

TEP signature of supra-threshold motor cortex activation
In the second part of the study, topographic analysis was used to compare TEPs elicited by sub-threshold and suprathreshold M1 stimulation with the aim of characterizing the TEP signature of supra-threshold activation of this cortical area.
We envisaged two sources of activity that might additionally contribute to the global topography of supra-threshold M1 TEPs. First, as supra-threshold stimulation of M1 produces a muscular contraction, these TEPs are probably interlaced with additional somatosensory-evoked activity triggered by the contracting hand muscles. Second, suprathreshold M1 stimulation possibly triggers additional brain activity in relation to the synchronous excitation of a large numbers of M1 pyramidal neurons. Although it is true that a small fraction of these neurons can fire even when a subthreshold stimulus is applied (Kujirai et al. 1993;Nakamura et al. 1997), we assumed that the corresponding brain response would be weak and therefore would have only a small impact on the topographic pattern of sub-threshold M1 TEPs.
In our dataset, the most remarkable difference between sub-threshold and supra-threshold M1 TEPs was the appearance of an additional early-latency N17 component following supra-threshold stimulation. This early response results in the stimulation of different muscle groups, thereby producing artifacts of specific distribution and magnitude. This might represent a confounding factor when comparing topographies of TEPs evoked by stimulation of different cortical targets, which we tried to avoid by carefully cleaning out the artifacts using ICA.
Late-latency TEP components. At later latencies, both M1 stimulation and AG stimulation triggered high-amplitude responses having very similar topographies, suggesting that they originate from similar cortical sources. This leads us to speculate that late TEP components N100 and P180 could predominantly reflect brain activity triggered by TMSinduced sensory stimulation, in agreement with the body of evidence accumulated in past years (Ahn and Fröhlich 2021;Freedberg et al. 2020;Rocchi et al. 2021;Ross et al. 2022). Importantly, the N100-P180 complex could be related at least in part to the so-called vertex potential which can be elicited by any transient sensory stimulus, regardless of sensory modality (Mouraux and Iannetti 2009). We also cannot exclude the possibility that cortical stimulation of a given region may generate cortical responses that are unspecific of the network targeted by TMS.
Such as the vertex potential, the P180 component of TEPs was symmetrically distributed over both hemispheres and maximal at Cz, suggesting its likely association with the sensory response. Such conclusion is further supported by results of a recent study that investigated the contribution of sensory evoked potentials to TEPs by exploring the association between the absolute stimulation intensity and TEP amplitude across all subjects (Niessen et al. 2021). They assumed that since stimulation intensity is defined relative to the individual rMT, the strength of direct cortical activation and by consequence the magnitude of corresponding TEPs should be similar across subjects and largely independent of the absolute TMS intensity. By contrast, stimulation of superficial structures may lead to sensory-evoked responses strongly depending on the absolute intensity. A significant regression between both variables was observed at latencies beyond 170 ms, leading the authors to attribute these late TEP components to the response to peripheral stimulation rather than to a direct cortical activation by TMS.
Unlike P180, and unlike the negative component of vertex potentials, the N100 TEP component identified in the present study was lateralized, indicating that additional cortical sources might be active at this time that are common for both TMS stimulation sites but not related to the sensory-evoked vertex potential. Our observations are in line with the findings of two previous studies using different experimental approaches. Rocchi et al. described a lateralized N100 component when comparing sub-threshold M1 TEPs with different sham conditions and various levels of sensory masking (Rocchi et al. 2021). However, they also both motor and sensory cortex in the stimulated hemisphere (Petrichella et al. 2017). Based on these observations, we hypothesize that the difference in P60 topography could be largely explained by the presence of a somatosensoryevoked potential triggered by the hand muscle contraction following supra-threshold stimulation.
As discussed earlier, the later TEP components N100 and P180 elicited by sub-threshold M1 stimulation can be mostly attributed to sensory-evoked vertex potential caused by the TMS discharge. If so, one would expect suprathreshold M1 stimuli to evoke a stronger response having the same topographic pattern. Yet, our data show clear topographic dissimilarity in the time intervals encompassing both components, with the t-map revealing additional negativity over the stimulated hemisphere in supra-threshold M1 TEPs. A similar observation was previously made by Gordon et al. (Gordon et al. 2018) who attributed this topographic difference to the somatosensory response to the evoked muscular contraction. Alternatively, we speculate that this late-latency activity specific to supra-threshold M1 stimulation might represent an inhibition of the ipsilateral sensorimotor cortex occurring in response to a strong, unexpected excitation of the M1. The N100 component of M1 TEPs has been previously linked to cortical motor inhibition as its amplitude was found increased when subject prepared to resist the TMS evoked contraction (Bonnard et al. 2009), decreased during movement (Kičić et al. 2008), decreased when contractions were actively facilitated (Nikulin et al. 2003) and correlated with the so-called cortical silent period (Farzan et al. 2013). It was also suggested that this component is at least partially mediated by slow metabotropic GABA B receptors, whose pharmacological activation can enhance the N100 amplitude (Premoli et al. 2014).

Usefulness and limitations of topographic analysis in TEP evaluation
In the present study, we show that the differential engagement of cortical sources contributing to TEPs can be successfully detected and characterized using topographic dissimilarity analysis, and that the spatiotemporal features of individual TEP components can be extracted and compared across datasets using microstate analysis. Together, our results allowed us to identify components of interest for future evaluation of the effect of experimental manipulation, thus demonstrating the practical utility of topographic analysis of TEPs.
The multiple advantages of topographic analysis have been already mentioned: it allows to consider all channels in a single measure, thus minimizing the loss of information; it identifies time intervals of interest in a purely datadriven manner, without the necessity of a prior choice; it has been described in other studies where it was referred to as the N15 or N18, and has been attributed to TMS activation of the ipsilateral M1 (Bonato et al. 2006) or the premotor cortex (Litvak et al. 2007). The peak-to-peak amplitude of the N15-P30 complex was found to correlate with the amplitude of MEPs (Mäki & Ilmoniemi, 2010). Therefore, this early component could be directly linked to the cortico-spinal output of M1. It is also possible that similar activity of much smaller magnitude occurs in response to sub-threshold M1 stimulation in a subset of participants, which would explain the period of topographic inconsistency that we identified in sub-threshold M1 TEPs at the latency corresponding to N17.
Almost all subsequent TEP components occurred at similar latencies for sub-threshold and supra-threshold stimulation. The P30 component had a greater amplitude for supra-threshold stimulation as compared to sub-threshold stimulation, but its mean topographical distribution was very similar for the two stimulation conditions, and the microstate analysis labelled it with the same microstate class. Therefore, the P30 could represent a cortical process that is triggered by both sub-and supra-threshold M1 stimulation. Its topography showed a clear positive peak at the proximity of the stimulation site and over the midline, suggesting activity originating from M1 and/or closely connected cortical areas, possibly including contralateral premotor or sensorimotor cortices, as suggested in previous studies (Esser, Huber et al. 2006;Litvak et al. 2007;Ahn and Fröhlich 2021). The P30 could thus reflect activity spread from M1 and could provide valuable information about the excitability of the M1 and its cortico-cortical connections.
Although peaking at similar latencies, almost all other TEP components elicited by sub-threshold and suprathreshold M1 stimulation exhibited markedly distinct topographies, indicating that they were generated by different combinations of cortical sources. It was previously established that motor re-afference may affect M1 TEPs starting from approximately 40 ms post-stimulus (Ilmoniemi and Kičić 2010;Paus et al. 2001). Coincidentally, the t-map of the P60 component revealed by TANOVA closely resembled the topography of the N20 component of somatosensoryevoked potentials elicited by electrical stimulation of the median nerve (Allison et al. 1989;Buchner et al. 1994;Van de Wassenberg et al. 2008). The latency of the two components would also match, considering that the earliest change in muscle spindle activity is typically observed around 35 ms after TMS stimulation of M1 (Rothwell et al. 1990). Finally, a previous study investigating M1 TEPs elicited by stimulation at the threshold intensity for eliciting a motor response reported a significantly larger P60 in TEPs accompanied by a MEP than in TEPs with no peripheral response. Moreover, the source of this component was localized in helps avoid the bias introduced by the arbitrary choice of the reference electrode. However, similarly to other methods employing point-by-point comparisons, topographic analysis assumes spatiotemporal consistency of TEP components across subjects -an assumption that has been disproved in the past (Ozdemir et al. 2021). While the topographic inter-subject variability is partially mitigated by the effect of volume conduction in EEG data, the variability in TEP latency might be problematic, particularly when considering the very early and very transient TEP components. Furthermore, since topographic analysis relies on normalized scalp voltage maps, it does not provide information on the strength of the observed response. It is therefore suitable to combine it with additional analyses to evaluate the amplitude of the elicited response.
At this point, we propose that the spatiotemporal features obtained by microstate analysis might also serve as a valuable prior for the extraction of the peak amplitude of individual TEP components at the subject level.

Conclusion
We have used topographic analysis to characterize the spatiotemporal features of TEPs and to identify TEP components most suitable to provide information about the properties of stimulated cortical areas and associated brain networks. By comparing TEPs from two different cortical areas, we confirmed that early TEP components until approximately 80 ms post stimulus likely represent, at least to a considerable extent, genuine TMS-evoked cortical activity. By contrast, later TEP components N100 and P180 showed site unspecific features, suggesting that they mostly reflect responses to the peripheral activation of sensory receptors, and as such should be interpreted with caution. Next, we evaluated the specific TEP signature of supra-threshold M1 activation as compared to sub-threshold M1 activation. We found that supra-threshold stimulation elicits an early N17 component which presumably corresponds to the activation of M1 neurons of the cortico-spinal tract. Topographical differences between sub-and supra-threshold M1 TEPs were also observed at later latencies, some of which may reflect somatosensory processing of inputs triggered by the muscular contraction. We conclude that topographic analysis might represent an advantageous method for TEP evaluation and interpretation. jurisdictional claims in published maps and institutional affiliations.
Springer Nature or its licensor holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.
Yaffe RB, Kerr MS, Damera S, Sarma SV, Inati SK, Zaghloul KA (2014) Reinstatement of distributed cortical oscillations occurs with precise spatiotemporal dynamics during successful memory retrieval. Proceedings of the National Academy of Sciences,111(52), 18727-18732 Publisher's Note Springer Nature remains neutral with regard to