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, Hallett, Rossini, Pascual-Leone, & Group, 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 pulses delivered manually using a MagPro X100 magnetic stimulator (MagPro X100 with MagOption; MagVenture, Farum, Denmark), inducing electrical current with anterior–posterior posterior–anterior (AP-PA) direction.
M1 stimulation. The stimulation target over the left M1 was defined as the most optimal position to elicit motor-evoked 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, Z'Graggen, Kohl, Rösler, & Kaelin-Lang, 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 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 & 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, Haufe, & Tangermann, 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.
It is necessary to emphasize that we took great care not to remove 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, peaks were identified in the average waveforms and named according to their polarity and approximate latency as measured at electrode Cz. Peak latencies were then extracted automatically from local maxima in the time course of the Global Field Power (GFP) that was calculated from the grand average waveform of each condition as the standard deviation of voltage across all channels and as a function of time (Equation 1). As such, GFP relays 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 & Skrandies, 1980).
Planned comparisons
To evaluate the effect of stimulation site, TEPs elicited by AG stimulation were compared to the TEPS elicited by sub-threshold 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 sub-threshold M1 stimulation should thus be more comparable to TEPs elicited by stimulation of non-motor cortical areas such as the 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 test-retest variability.
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, Kottlow, Stein, & Melie-García, 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, Weusmann, Stein, and Koenig (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; sub-threshold vs supra-threshold M1 stimulation), we used the measure of Global Dissimilarity (DISS, Equation 2) also introduced by (Lehmann & 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 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 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. First, the optimal number n of microstate classes was determined using the cross-validation method proposed by (Koenig, Stein, Grieder, & Kottlow, 2014). Individual data were randomly split in half into a training and a testing dataset, and averaged. Then, the training TEP waveforms were used to compute microstate models with the number of classes ranging from 3 to 12. The mean spatial correlation of the testing dataset with each microstate model was computed according to Equation 3. This process was repeated 5000 times for variable subsets of the data. The mean correlation from all runs was averaged and compared across models to indicate the one yielding the maximum mean correlation, and the corresponding number n of microstate classes was retained for the subsequent analysis.
In the next step, n class maps were identified in a dataset created by concatenating the group-level average TEP waveforms of all compared stimulation conditions (AG stimulation, sub-threshold 1 stimulation and supra-threshold M1 stimulation). 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 TEP waveforms into non-overlapping intervals of dominance of one of the n class maps. For this purpose, we chose to use the k-means clustering algorithm (Murray et al., 2008) with 5000 permutations. In each cycle, the algorithm draws n random samples (topographies at n different timepoints) from the concatenated dataset and computes the spatial correlation between each sample and all the remaining topographies in the dataset. In the next step, the Global Explained Variance (GEV, Equation 4) was computed for each class map based on the correlation values. All class maps were subsequently re-defined by averaging topographies over all the timepoints where each given class map yielded the highest GEV as compared to other class maps. The spatial correlation and GEV were re-calculated and the whole process was repeated until no increase in GEV is be observed. The maximum obtained GEV was then associated with the final version of the obtained class maps, and the next iteration began with yet another random sample of n topographies. Finally, the set of n class maps that yielded the highest overall GEV 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.