A Structured ICA-based Process for Removing Auditory Evoked Potentials Reveals TMS-evoked Potentials and TMS-modulated Oscillations

Jessica M. Ross (  jross4@bidmc.harvard.edu ) Harvard Medical School Recep A. Ozdemir Harvard Medical School Shu Jing Lian Berenson-Allen Center for Noninvasive Brain Stimulation, Beth Israel Deaconess Medical Center Peter J. Fried Harvard Medical School Eva M. Schmitt Hinda and Arthur Marcus Institute for Aging Research, and Deanna and Sidney Wolk Center for Memory Health, Hebrew SeniorLife, Boston Sharon K. Inouye Harvard Medical School Alvaro Pascual-Leone Harvard Medical School Mouhsin M. Sha Harvard Medical School


Introduction
Transcranial magnetic stimulation (TMS) has gained increasing importance and application in the neurophysiologic characterization of healthy aging [1][2] and a variety of neurological disorders, including Alzheimer's disease [1][2], epilepsy [3], autism spectrum disorders [2], and schizophrenia (see [4] for a we used non-M1 stimulation sites for the active stimulation conditions. We quantify the changes that occur when removing AEP using this method and characterize the components that account for TEP nonspeci city. This manuscript presents these four speci c contributions: 1) We assess the utility of this approach in two populations: younger and older adults, collected with differing sham procedures and TMS-EEG systems. 2) We apply a group level analysis of the impact of removing AEP with this method. 3) We quantify the speci city of the TEPs across different non-motor stimulation sites and between individuals before and after this procedure. Finally, 4) we describe the residual activity. Our hypothesis was that effectively removing AEP would reveal stimulation site-speci c and participant-speci c responses in early, mid latency and late time windows.

Participants
Data used in the present analysis were from two ongoing TMS-EEG studies at the Berenson-Allen Center for Noninvasive Brain Stimulation at Beth Israel Deaconess Medical Center. All participants provided written informed consent before enrollment according to the Declaration of Helsinki. Both studies were approved by the Institutional Review Board of the Beth Israel Deaconess Medical Center, Boston, MA. The rst cohort consists of 10 healthy adult control participants (9M/1F, age = 42.2 ± 18.8 yrs., range = 19 to 65) from a TMS-EEG study of epilepsy [3]. The second cohort consisted of 24 older adults, collected as part of a study of postoperative delirium called The Successful Aging after Elective Surgery study (SAGES; 10M/14F, age = 72.0 ± 6.6 yrs., range = 65 to 89) [38]. Participants selected for SAGES did not have clinical dementia, were scheduled for major surgery, and we used pre-surgery baseline TMS-EEG visit data for this analysis. We replicated and validated the ndings from the rst cohort with the second cohort. The two cohorts represent different age groups, sham TMS protocols, and TMS-EEG systems, and are therefore ideal for addressing the robustness and strengths/weaknesses of this structured ICA-based approached to isolating AEP from TEP independent of experimental design.

Transcranial magnetic brain stimulation set-up and procedure
In the younger cohort, TMS was delivered using a Nexstim eXimia stimulator with real-time MRI neuronavigation (NBS software v3. 2.1;Nexstim,Helsinki,Finland), following standard protocols [39]. Single TMS pulses were applied to left dorsolateral prefrontal cortex (L DLPFC), left intra-parietal lobule (L IPL) and a sham condition to left M1 delivered with an active coil tilted away from the scalp at 90 o .
DLPFC and IPL were anatomically de ned, respectively, as the superior half of the middle frontal gyrus approximately 3 cm anterior to precentral sulcus and the superior edge of angular gyrus roughly 1 cm inferior to intraparietal sulcus. Coil orientation was set with coil handle perpendicular to the targeted gyri inducing a biphasic anterior-posterior-posterior-to-anterior current in the underlying cortex. No auditory noise masking or electrical stimulation was used. Participants were asked to wear earplugs during sham and active-TMS trials to protect their hearing. Throughout the visit, participants were seated in an adjustable chair. The motor hotspot for eliciting motor evoked potentials (MEPs) in the right rst dorsal interosseous (FDI) muscle was determined by delivering single pulse TMS to the hand knob of left motor cortex with the coil tangential to the scalp and oriented orthogonal to the motor strip (~45 degrees from anterior-posterior midline), and searching around until nding the site where the largest and most consistent MEPs were produced in FDI. Resting motor threshold (RMT) was then determined as the minimum stimulation intensity (% of maximum stimulator output) needed to elicit MEPs ≥ 50 µV in the relaxed FDI in at least 5 out of 10 pulses. Following determination of RMT, 120 single pulses of TMS were applied at 120% RMT with randomly jittered inter-stimulus intervals (3-5 seconds). TMS operators monitored participants for wakefulness. See Vink et al. [40] for more details.
In the older SAGES cohort, TMS was delivered using a gure-of-eight shaped coil with dynamic uid cooling, a biphasic waveform, and anterior-posterior-posterior-anterior current direction in the brain (MagVenture Cool-B65) attached to a Magpro X100 stimulator, following standard protocols [39].
MagVenture Cool-B65 A/P coil was used with a 2.7 cm thick plastic spacer for sham stimulation and white noise masking was presented through earplug-earbuds at the maximum volume comfortable for each participant. Although auditory noise masking was used in the older cohort, 20 participants reported hearing the TMS "click" and had AEPs. These 20 participants were used for all pre/post AEP removal comparisons (7M/13F, age = 71.9 ± 5.8 yrs., range = 65 to 86). Small current pulses (between 2-4 mA) were delivered over the left forehead, over the frontalis muscle, using surface electrodes (Ambu Neuroline 715 12/Pouch) to approximate somatosensory sensations arising from skin mechanoreceptors and scalp muscles during the active-TMS condition. The motor hotspot for eliciting motor evoked potentials (MEPs) in the right FDI muscle was determined by methodically delivering single pulse TMS to targets of the hand knob of left motor cortex with the coil tangential to the scalp and oriented orthogonal to the central sulcus/motor strip (~45 degrees from mid-sagittal plane) until the largest and most consistent MEPs were produced in FDI. RMT was determined as the minimum stimulation intensity needed to elicit MEPs ≥ 50 µV in the relaxed FDI in at least 5 out of 10 trials. Following determination of the FDI hotspot and RMT, 120 single pulses of TMS were applied at 120% RMT with randomly jittered inter-stimulus intervals (3-5 seconds) to the DLPFC and IPL targets, and sham TMS over M1. To select the non-motor targets in the SAGES cohort, group-level resting-state functional networks maps were used based on a 7 networks parcellation covering cortical and subcortical structures [41]. Con dence maps for each RSN were used, representing the con dence of each vertex belonging to its assigned network across a sample of 1000 healthy participants (expressed as valued between -1 and 1), with larger values indicating higher con dence. By using group-level functional parcellations and con dence maps on the Montreal Neurological Institute (MNI) template brain, we were able to target the most consistent and reliable regions within the left angular gyrus (L-IPL: - 55.1, -70.5, 27.7) and dorsolateral prefrontal cortex (L-DLPFC: -50.5, 30.4, 31.8) that had the highest likelihood of occurring in the default mode and frontoparietal resting-state networks, respectively. A custom processing pipeline was developed to take each participant's anatomical MRI, create a non-linear transform from the participant's native space to MNI space and then use the inverse of that transform to bring back the coordinates into participant's space using FSL's FNIRT tool. The transformed coordinates along with individual high-resolution T1w images were then imported into the Brainsight™ TMS Frameless Navigation system (Rogue Research Inc., Montreal, Canada), and co-registered to digitized anatomical landmarks for online monitoring of coil positioning. TMS operators monitored participants for wakefulness.

Electroencephalography acquisition and preprocessing
In the younger cohort, whole scalp 60-channel EEG was collected with a TMS-compatible ampli er (eXimia EEG, Nexstim Ltd) and labeled according to the extended 10-20 international system. EEG data were online referenced to an additional electrode on the forehead. Impedances were maintained below 5kΩ at a sampling rate of 1450 Hz. EEG signals were digitized using a Nexstim ampli er. Digitized EEG electrode locations on the scalp were also co-registered to individual MRI scans using Nexstim Brainsight™ TMS Frameless Navigation system. See Vink et al. [40] for more details.
In the SAGES cohort, whole scalp 64-channel EEG was collected with a TMS-compatible ampli er (actiCHamp, Brain Products GmbH, Munich, Germany) and labeled according to the extended 10-20 international system. EEG data were online referenced to AFz. Impedances were maintained below 5kΩ at a sampling rate of 5000 Hz. EEG signals were digitized using a BrainVision ampli er (BrainCHamp,BrainVision Recorder software,v. 1.21). Digitized EEG electrode locations on the scalp were co-registered to individual MRI scans using Nexstim Brainsight™ TMS Frameless Navigation system. EEG data were pre-processed o ine using custom scripts and EEGLab v2019.0 [42] in MATLAB R2019a (Mathworks, Natick, MA, USA). Data were segmented into 3000 ms epochs, each starting 1000 ms before (pre-stimulus) and ending 2000 ms (post-stimulus) following TMS pulse, respectively. Baseline correction was performed by subtracting the mean pre-stimulus (-900 to -100) signal amplitude from the rest of the epoch in each channel. Following baseline correction, data were visually inspected to identify and remove extremely noisy channels with high amplitudes. Zero-padding between -2 ms and 14 ms time range were then applied to remove the early TMS pulse artifact from the EEG data. All zero padded epochs were then tagged based on voltage (≥100 μV), kurtosis (≥ 3), and joint probability (single channel-based threshold ≥ 3.5sd; all channel-based threshold ≥ 5sd) metrics to identify excessively noisy epochs. Visual inspection was performed on the tagged epochs for the nal decision on removal of noisy epochs.
An initial round of fast independent component analysis (fICA) [43] was performed to identify and remove exclusively components with early TMS evoked high amplitude electrode and EMG artifacts [35,44]. After the rst round of fICA, the EEG data were interpolated for previously zero-padded time window around TMS pulse (-2 ms to 14 ms) using linear interpolation, band pass ltered using a forwardbackward 4th order Butterworth lter from 1 to 100 Hz, notch ltered between 57 and 63 Hz, and referenced to global average. Missing/removed channels were interpolated using spherical interpolation, and DLPFC, IPL and sham stimulation blocks were merged, in that order. Because recordings were made with 60/64 channels, and the signals were unlikely to have that many independent sources, PCA was used to reduce dimensionality prior to ICA. This approach can improve decomposition [45][46] and signal to noise ratio of large sources [47]. Data were reduced into no fewer than 30 dimensions in the younger cohort (median=32. 5, range=30-50), and no fewer than 35 dimensions in the older cohort (median=45, range=35-60), optimized at the individual level to maximize uniqueness and dipolarity of resulting components while minimizing residual noise components [47]. A second round of fICA was run on the merged data to manually remove all remaining artifact components [44] including eye movement/blink, muscle noise (EMG), single electrode noise, TMS evoked muscle, and cardiac beats (EKG). The data were low pass ltered with a 4th order Butterworth lter at 50 Hz.
Components that met the following criteria were classi ed as AEP components and removed: (1) the time-course has three peaks, P50-N100-P200, with the lowest amplitude in the P50 compared with N100 and P200; (2) the scalp topography re ects left/right symmetry and a central distribution anterior to Cz; and (3) the component is shared across both active stimulation sites and sham stimulation. See Discussion for more details about AEP classi cation criteria. Conditions were then unmerged for subsequent analyses.
In both rounds of fICA, a semi-automated artifact detection algorithm incorporated into the open source TMS-EEG Signal Analyzer (TESA v1.1.0-beta) extension for EEGLAB was used to classify and visually inspect components based on their frequency, activity power spectrum, amplitude, scalp topography, and time course [48] (http://nigelrogasch.github.io/TESA/).

Analysis of TMS-evoked potentials
Time-course and topography: Latency and amplitude of peaks in the pre-and post-AEP removal TMSevoked cortical response were identi ed using an automated peak nder algorithm [49] and scalp topographies were plotted at those peak latencies.

GMFA/LMFA (and windows):
Global mean eld potential (GMFA) across all channels was used to quantify global synchrony and excitability [50][51] in the 14-400 ms post-TMS window. We computed Global Mean Field Power (GMFP) following: Where V i (t) is the voltage at electrode i at a certain point in time, V mean(t) is the mean of instantaneous TEP across electrodes, and K is the number of electrodes.
Local mean eld potentials (LMFA) from electrodes near the active stimulation sites were used to quantify local activity under the coil [50,[52][53] in the 14-400 ms post-TMS window. Channels included in each region of interest were determined as all electrodes bordering stimulation target locations. For the younger cohort, the DLPFC ROI included AF1, F7, F5, F1, FT7, FC5, and FC3 and the IPL ROI included TP7, CP5, CP3, P7, P3, and PO3. In the older cohort, the DLPFC ROI included AF3, F7, F5, F3, FT7, FC5, and FC3, and the IPL ROI included TP7, CP5, CP3, P7, P5, P3, PO7, and PO3. See Supplementary Figure S1 for electrode arrays color coded to show LMFA ROIs for younger adult (S1A) and older SAGES (S1B) cohorts. In addition, GMFA and LMFA were calculated for three windows corresponding to the three subcomponents of AEP. Windows were de ned in each study cohort for GMFA, LMFA of the DLPFC ROI, and LMFA of the IPL ROI, using the latency of the minimum amplitude between peaks in the sham stimulation condition (determined using the peak nder() function, which nds peaks with amplitudes that are greater than two standard deviations above baseline [49]). See Figure 2D for time-windows de ned for GMFA analysis in the younger cohort, and Supplementary Figure  Similarity Index: Similarity between TEPs across site and between subjects at each site, pre and post AEP removal, was quanti ed in vector space using the cosine similarity index [12,20]. We rst generated a TEP matrix for each subject (from averaged responses) with a xed window size (385ms) covering EEG responses from 14 to 400ms. Each TEP matrix contains millisecond voltage values from all channels with a 63x385 matrix size. We then used cosine similarity to quantify similarity index (SI) between matrices following: Where SI XY is the cosine similarity between TEP matrices x and y for a given stimulation site, and n is the number of channels, X it and Y it are the i th vector of all channels at time t. Within-subject site SI was computed as the mean SI from the diagonal of site comparison grids. Between-subject SI was computed as the mean SI of half grid, excluding the diagonal, of within-site grids (DLPFC by DLPFC, etc.). To enable assessment of similarity between the DLPFC and IPL datasets, SI was also calculated for three timewindows de ned by residual TEP in active stimulation conditions post-AEP removal: early, mid, and late latency. These windows were de ned using GMFA area under the curve post-AEP removal active stimulation conditions averaged across participants.
We performed multiple univariate tests to describe how the AEP-removal method affects each outcome variable separately. Multiple univariate tests are preferred over one multivariate MANOVA here for a number of reasons [54]. First, our research question dictates a model that tests for effects of the method on each outcome variable independently, with no particular interest in a linear composite of the outcome variables. Our approach best describes how using this AEP-removal method will practically in uence a variety of datasets if applied. We acknowledge that redundant information is being obtained using multiple ANOVAs but suggest that this redundancy supports that removing AEP using this method has robust effects on mean eld potentials across ROIs and, particularly at mid and later latencies. However, to address concerns about experiment-wise Type 1 error for performing multiple ANOVAs, we additionally use an additive Bonferroni inequality. For m = 12 tests, α of 0.05 becomes 0.004 (0.05/12).

EEG Source Reconstruction:
All TMS evoked EEG source reconstruction was performed using Brainstorm [55]. First, digitized EEG channel locations and anatomical landmarks of each subject were extracted from Brainsight™ (nasion 'NAS', left pre-auricular 'LPA', and right pre-auricular 'RPA' points), and registered onto individual MRI scans in brainstorm. Next, the EEG epochs, -500ms to 1000ms with respect to TMS pulse for each TMS trial were uploaded, and the average epoch time series was generated for each subject. Forward modeling of neuro electric elds was performed using the open MEEG symmetric boundary element method, all with default parameter settings. Noise covariance was estimated from individual trials using the pre-TMS (-500 to -1ms) time window as baseline. The inverse modeling of the cortical sources was performed using the minimum norm estimation (MNE) method with dynamic statistical parametric mapping (dSPM) and constraining source dipoles to the cortical surface. The resulting output of EEG source reconstruction was the MNE current density time series for each cortical vertex.
Residual Mid-Latency Components: As a follow-up we examined the residual TEP in a mid-latency window. We used two approaches in an attempt to isolate the mid-latency components in sham-(1) a percent variance threshold to identify components contributing the most variance in this mid-latency window, and (2) non-parametric permutation based suprathreshold cluster size analysis [56]. We looked for clusters in the post pulse GMFP that were greater than the baseline period, and we determined which components contributed the highest percent variance to each of these clusters. See Results 3.4 for more details. Due to redundancy in the results of these two approaches, only the percent variance threshold approach (1) is described here. See Supplementary material for the results of the non-parametric permutation-based cluster size analysis (2).

Results
Results of all analyses are shown for the younger cohort in and summarized results are shown for the older cohort in Figure 6. See Supplementary for additional Tables and Figures for both  cohorts. 3.1 The AEP All 10 participants in the younger cohort had an AEP that met all three criteria (time-course, topography, and shared across conditions including sham). See Figure 1A for a representative participant's TEP and scalp topography of the component classi ed as AEP in this participant's data. The ERP image of this component from the merged conditions shows a shared time-course and polarity of AEP across active and sham stimulation conditions ( Figure 1B for same representative participant). Averaged TEP and topography plots across all AEP components from the full cohort show the expected stereotypical timecourse and topography ( Figure 1C). At the group level, AEP components (14-400ms) showed low withinsubject site-speci city with an averaged between-site SI of 0.83 ± 0.04. The AEP components showed moderate subject speci city, with an averaged within-site between-subject SI of 0.49 ± 0.03. For all between-and within-site SIs, see Supplementary Table S1.
In the older cohort, 20 of 24 participants had an AEP that met all three criteria. AEP components (14-400ms) again had high similarity across sites, with an averaged within-subject between-site SI of 0.66 ± 0.06. Between-subject similarity was again moderate, with an averaged within-site between-subject SI of 0.36 ± 0.02. For all between-and within-site SIs, see Supplementary Table S1. See Figure 6A for the averaged AEP component time-course and topography from active stimulation of the DLPFC target in the older cohort.

Effects on the TEP
Time-course and topography: Pre/post AEP removal comparisons at the group level revealed an amplitude reduction that was most noticeable at approximately 100 and 200 ms latencies. See Figure 2A-C for TEPs and scalp topographies before and after AEP removal in the rst cohort (N=10). TEP from stimulation of DLPFC ( Figure 2A) had peaks at 28, 44,63,110,188,308, and 363 ms pre-AEP removal, and AEP removal revealed smaller peaks at 123 and 152 ms. See Figure 2B IPL stimulation at 15,47,76,108,188,304, and 360 ms pre-AEP removal and the revealed peak at 117 ms post-AEP removal, and Figure 2C sham stimulation at 31,99,188,306, and 370 ms pre-AEP removal and the revealed peaks at 57 and 177 ms post-AEP removal. Scalp topographies in this mid-latency window shifted to posterior distributions with removal of AEP.
Earlier peaks appeared to be stimulation site speci c and have unique topographies, less effected by removal of AEP than later peaks.
In the older cohort, with removal of AEP there was a similar amplitude reduction in mid-latency peaks of the TEPs and revealing of smaller peaks in that window with a posterior shift in topographies. See Figure  6B for pre-and post-AEP removal TEPs and scalp topographies in the SAGES cohort with DLPFC stimulation (note pre/post changes at 117 and 193 ms).
In the older cohort, in the full 14-400 ms time window, we found that removing AEP reduces GMFA AUC (F (1,19)  Using the peak/trough analysis, we found windows de ned for the DLPFC ROI as 14 to 54 ms, 54 to 146 ms, and 146 to 400 ms, and for the IPL ROI as 14 to 52 ms, 52 to 144 ms, and 144 to 400 ms (Supplementary Figure S2B-C). Removing AEP did not reduce LMFA area under the curve in the P50 window in either ROI (DLPFC ROI: F (1,9) = 0.47,p = 0.51;IPL ROI: F(1,9)  To address concerns about experiment-wise Type 1 error for performing multiple ANOVAs, we additionally use an additive Bonferroni inequality. As described in Methods, for m = 12 tests, α of 0.05 becomes 0.004 (0.05/12). Our general pattern of results remains signi cant across measures and cohorts, showing a reduction in GMFA and LMFA AUC for full TEP, and with most robust changes in the N100 and P200 latency windows. Interactions between version and site do not remain signi cant for most outcome measures, excepting GMFA of N100 in the older cohort.
See Supplementary Table S2 for results of all GMFA and LMFA analyses, with test statistics, using α of 0.01, 0.05, and 0.004.

Similarity Index:
In the younger cohort, in full TEPs (14-400 ms), we found a reduction in between-subject SI pre-(0.48 ± 0.012) to post-(0.15 ± 0.013) AEP removal (t(2) = 17.82, p < 0.0001), indicating increasing subject speci city within stimulation site with removal of AEP. We found a reduction in between-site SI pre-(0.71 ± 0.019) to post-(0.47 ± 0.018) AEP removal (t(2) = 9.17, p < 0.0001), indicating an increase in site speci city with removal of AEP, although it should be noted that 0.47 ± 0.018 was still moderately similar across site. See Figure 3A-B for between-subject pre and post AEP removal grids (A) and between-subject and between-site SI averages pre and post AEP removal (B).
In the active stimulation conditions, GMFA post AEP removal appeared to have three latency windows, separated by troughs at 77.24 and 273.79 ms in DLPFC and 77.24 and 270.34 ms in IPL, with trough latencies averaged across active sites at 77.24 and 271.72 ms. These windows were chosen to enable assessment of similarity in residual TEP (post AEP removal) and were de ned using peak nder() [49] as early: 14 to 77.24 ms, mid: 77.24 to 271.72 ms, and late: 271.72 to 400 ms ( Figure 3C). In the early window, between-site SI was 0.39 ± 0.076 pre-AEP removal, indicating site speci city, and 0.24 ± 0.056 post-AEP removal with no signi cant change from pre-to post-AEP removal (t(2) = 1.61, p = 0.092). In the mid-latency window, between-site SI was 0.81 ± 0.023 pre-AEP removal, indicating high similarity across sites in this window, and reduced to 0.57 ± 0.036 post-AEP removal (t(2) = 5.68, p = 0.0024). In the late window, between-site SI was 0.42 ± 0.05 pre-AEP removal and 0.35 ± 0.034 post-AEP removal with no signi cant change from pre-to post-AEP removal (t(2) = 1.15, p = 0.16). See Figure 3D for between-site SI averages pre and post AEP removal in the three latency windows.

See Supplementary Table S3 for mean between-subject and within-subject between-site SIs for all stimulation conditions.
Similarly, in the older cohort, in full TEPs (14-400 ms), we found a reduction in between-subject SI pre-(0.37 ± 0.036) to post-(0.16 ± 0.019) AEP removal (t(2) = 5.03, p = 0.0037), and a reduction in betweensite SI pre-(0.58 ± 0.027) to post-(0.37 ± 0.028) AEP removal (t(2) = 5.29, p = 0.0031). As described above, AEP components (14-400ms) had an averaged between-subject SI of 0.36 ± 0.02 and an averaged between-site SI of 0.66 ± 0.06-the almost identical SIs of AEP components and of pre-AEP removal TEP (but reduced in post-AEP removal TEP) indicates that the AEP components are dominating the signal-tonoise ratio of the evoked response prior to removal.
In the active stimulation conditions, GMFA post AEP removal showed three latency windows, separated by trough latencies averaged across active sites at 86.80 and 285.20 ms. The three latency windows were de ned using peak nder() [49] as early: 14 to 86.80 ms, mid: 86.80 to 285.20 ms, and late: 285.20 to 400 ms. In the early window, between-site SI was 0.27 ± 0.021 pre-AEP removal and remained low, 0.19 ± 0.038, post-AEP removal (t(2) = 1.98, p = 0.060). In the mid-latency window, between-site SI reduced from 0.72 ± 0.030 to 0.51 ± 0.040 pre-to post-AEP removal (t(2) = 4.17, p = 0.0070). In the late window, between-site SI was 0.24 ± 0.017 pre-AEP removal and 0.21 ± 0.020 with no change pre-to post-AEP removal (t(2) = 1.00, p = 0.19). See Figure 6C for between-subject and between-site SI averages pre and post AEP removal in the older cohort, and Supplementary Table S3 for mean between-subject and withinsubject between-site SIs for all stimulation conditions.

Source Analysis:
When AEP is present, topographical distribution of selected peaks in TEPs in both datasets have similar spatial characteristics compatible with the signal morphology and topography of AEPs reported in the recent TMS-EEG and EEG literature [35,[57][58][59]. While the presence of AEPs resulted in potentials with a uniform spatial topography and source activations at the expected time points, site-speci city of the evoked response and distinct source activations became clearly visible after removing AEP. See Figure 4 for spatial and temporal characteristics of TEPs pre-and post-AEP removal for a representative subject.

Residual Mid Latency Components
Removing AEP enables further exploration of residual TEP, so as a follow-up we examined the residual TEP that exhibits moderately low site speci city in the mid latency window. Although between-site similarity was reduced in this mid latency window, it was still higher than the midpoint of SI range (0 to 1), at 0.57 ± 0.036 in the younger cohort and 0.51 ± 0.040 in the older cohort, indicating some residual similarity across stimulation sites. One possibility is that the AEP was not being fully removed. However, the altered TEP time course and posterior shift in scalp topography in this mid latency window (Figure 2A-C and 6B) suggest that removing AEP revealed smaller amplitude components that are unique from AEP.
In an attempt to isolate the mid latency components in sham, we used a percent variance threshold to identify components contributing the most variance in this mid latency window. We looked for clusters in the post pulse GMFP that were greater than the baseline period, and we determined which components contributed the highest percent variance to each of these clusters. Percent variance threshold for components that contributed most to the mid latency window: Percent variance (PVAF) was calculated for the mid latency window for each participant. See Figure 5A for PVAFs of all components contributing to the mid latency window in one representative participant. In most participants, we found 1-2 components with a peak PVAF above 28%. This 28% threshold was selected using a data driven approach to isolating sources contributing most to the mid latency window.
The representative participant shown in Figure 5 had a peak PVAF in component 2. The scalp topography of IC2 (Fig. 5B) had a posterior distribution characteristic of occipital alpha. The time-course of IC2 in all trials (Fig. 5C) showed presence in baseline, increase in amplitude after the TMS pulse, and oscillations at approximately 10 Hz (Fig. 5C-D). After removing IC2, GMFA of the sham condition, calculated across all channels, decreased in both the pre-and post-stimulation windows (Fig. 5E). All peak PVAF components across participants that exceed 28% variance (N=7) showed activity in the pre-TMS baseline window, highest power oscillation in the alpha band, and GMFA reduction when removed in both the preand post-TMS windows. See Supplementary gure S3A for individual subject component topographies that exceeded 28% in this mid latency window. All these components across participants had a posterior distribution in the scalp topography, although some were at the midline and some suggested left/right laterality. See Figure 5F for the average TEP and topographies of all peak PVAF components that met the 28% threshold. Time-courses and topographies of the mid-latency components in sham suggested brain components that did not meet our time course/topography criteria for AEP. The components appeared to be occipital alpha with modulated power in the post-TMS window.
Between-subject and within-subject site similarity SI was calculated for the identi ed mid-latency components. Within-subject between-site SI was 0.53 ± 0.039, and between-subject SI was 0.12 ± 0.034 for these components. Together, these SI values indicated that the mid latency components did not exhibit site speci city but did exhibit subject speci city. Time course of component amplitude peaks showed no consistency across site within-subject, suggesting these components were not time locked to the TMS pulse.
Similarly, in the older cohort, we identi ed peak PVAF contributors in mid latency sham using the same 28% variance threshold. 16 participants had at least one component that met this threshold. See Figure  6D for the average TEP and topographies. These components exhibited posterior, either central or left/right lateral, scalp distributions, activity in the baseline period, highest power in the alpha band, and GMFA reduction in both the pre-and post-TMS windows when removed. See Supplementary Figure S3B for individual subject component topographies that exceeded 28% in this mid latency window.
Within-subject between-site SI was 0.32 ± 0.019, and between-subject SI was very low for these components, 0.10 ± 0.026. These SI values indicated that the mid-latency components did not exhibit high site speci city but did exhibit subject speci city. Time-course of component amplitude peaks between sites suggested these components were not time locked to the TMS pulse.
See Supplementary Table S4 for similarity of mid latency components between-subject and withinsubject between-site from all stimulation conditions, and Supplementary Table S5 for SI of TEPs pre-and post-removal of mid-latency components.

Discussion
In the present study, we show that auditory-evoked potentials (AEPs) evoked by the sound associated with the discharge of each TMS pulse can be isolated and extracted from the TEP after performing ICA on merged active and sham stimulation conditions. Using data from two separate studies, we show effective and conservative extraction of AEP with younger adult as well as an older population, with variations of sham stimulation protocols, and with different TMS devices and EEG systems. We show this method preserves residual TEP in early, mid and late latency windows that can be subsequently evaluated for speci city of stimulation site, group, and individual. The method not only preserves transcranial potentials but reveals non-speci c alpha modulations that were previously obscured by AEP. Our group level analyses of pre-and post-AEP removal suggest that the TEP is composed of transcranial evoked potentials, sensory-evoked potentials, and site-independent modulation of ongoing oscillations.
One aspect of our study that is unique is that the two cohorts that were used differed in a number of factors, including age, gender, the details of the sham application, and the use of auditory noise masking.
The consistency of results across the two cohorts demonstrates the robustness of the AEP removal method to these multiple experimental factors. E cacy and advantages of an ICA-based approach to removing AEP Options are needed for removing AEP from TEP to reveal other brain responses to TMS. Although earplugs and masking [6,17,31] can be used to attenuate AEP, and foam padding can reduce bone conduction of the sound [32], these techniques do not always work as effectively as anticipated. Many groups have observed AEP even after using these techniques [15][16][32][33][34]. Consistent with this literature, in our older cohort, in which auditory noise masking was used, a majority of participants perceived the 'click' and had an AEP. This may be due to contrasting acoustic properties of the 'click' sound and noise masking. However, regardless of the reason for the persistent AEP, options for effective handling of the AEP after data collection could be of great use, and we present one option that is effective. Further, we show that this method is also conservative in that it preserves the earliest TEP that is stimulation site speci c and preserves and reveals later activity that could carry information about subject-speci c neural modulations.
We show that removing AEPs using this structured ICA-based process signi cantly reduced GMFA (whole scalp) and LMFA (of left DLPFC and IPL ROIs) in the post-stimulation TEP (14 to 400 ms), driven by time windows consistent with the N100 and P200 temporal characteristics of AEP. Further, supported by cosine similarity analysis, we show that removing AEPs reduces TEP similarity between-subjects and between stimulation conditions. Similarity is reduced most in a mid-latency window, but nevertheless remains higher than mid-range for the 0-1 SI rating. Residual TEP in this window has a time course and topography unique from AEP, and follow-up analyses suggest this is alpha modulation that is not stimulation site speci c but is unique to individual subject.
One property of our isolated AEP components that might seem counterintuitive is that there are slight differences in amplitude of the AEP between sites ( Figure 1C). The acoustic properties of the 'click' were identical across conditions [6,17], but the AEP amplitudes were not. This could be due to differences between conditions in coil distance from the ears and/or bone conduction. Nikouline et al. [22] observed an amplitude difference between AEPs when the coil was pressed against the scalp, held 2 cm above the scalp in the air, or 2 cm above the scalp but with the use of a plastic spacer. In that study, all stimulation conditions were over the same left hemisphere M1 target, and time course of the AEPs were comparable.
The amplitude of N100-and P180 peaks were smallest when the coil was held above the scalp, intermediate when a spacer was used, and largest when the coil was pressed directly to the scalp, demonstrating how coil distance from scalp and bone conduction can contribute greatly to AEP amplitudes. Due to these results, AEP should not be expected to have comparable peak amplitudes under different stimulation conditions. Consequently, merely subtracting out the mean sham-evoked potential from active stimulation TEPs would not su ce in accounting for differential AEP amplitudes with other stimulation sites. In contrast, one of the strengths of the ICA approach is that as the underlying neural generators (and thus the resulting scalp EEG topographies) are stable, differential activations on distinct trials can be isolated and removed.

Classi cation of AEP Components
Features of the AEP time course, topography and presence in sham can be used for identifying and classifying AEP components, based on the literature from TMS-EEG and from other studies of auditory perception [60]. Because the auditory system is particularly sensitive to timing, this is re ected in temporal consistency of neural responses to sound. Auditory activity can be observed through multiple stages of the ascending auditory pathway with temporal precision, beginning with the auditory brainstem response (ABR; Jewett waves I-VII; 0-10 ms after stimulus onset) and followed by MLRs originating from medial geniculate nucleus and primary auditory cortex ("Mid-latency responses"; 10-60 ms after stimulus onset) [61]. Sensory AEP is generally described 40-200 ms after the stimulus, before P300 and other cognitive components [22,[26][27][62][63]. Because of axonal divergence up the auditory pathway, neuronal population increases and so do component amplitudes, with the N100 and P200 having the largest peak amplitudes. The AEP consists of at least three subcomponents, described as the "P50-N100-P160 complex" [64-66] and the "P1-N1-P2 complex" [67][68]. We use nomenclature consistent with TMS-EEG literature [18,25], and describe the subcomponents at 50, 100, and 200 ms of the well-documented sensory AEP. Because of the consistent time course of these subcomponents, with a smaller amplitude peak at 50 ms and larger amplitude peaks at 100 and 200 ms, time course can be used to identify AEP [27]. Our criteria for AEP classi cation adhered to this time course of three subcomponents at 50, 100 and 200 ms.
Source analyses of AEP subcomponents reveal that P50 originates in primary auditory cortex and N100 and P200 subcomponents originate in surrounding belt areas of A1 [27,[69][70]. Because of the mirrored, time-locked bilateral sources (left/right auditory cortices), distribution in scalp electrodes is symmetrical and gives the appearance of a single central deep dipole [29][30]. Our criteria for AEP classi cation included this expected central and symmetrical scalp topography.
Sham TMS protocols that use an active coil should result in a similar or identical AEP response if the acoustic properties of the "click" are the same across sham and active protocols [35]. Assuming the sham protocol is not itself inducing activation of cortex, the AEP response in sham represents non-TMSevoked potentials. Our criteria for selecting AEP components included a rigid adherence to this assumption-we only classi ed components as AEP if they were shared across stimulation site and sham stimulation.
Site-speci city of TMS-evoked potentials and other modulated oscillations Although there is evidence to support that TEP can have stimulation site speci city [9,[20][21][22], there is accumulating evidence for non-speci city in mid latency and late latency windows [15][16]20,24]. The present data provide support for the notion that the early responses to TMS are highly speci c to the site of stimulation, and thus likely represent transcranial-evoked elements. In contrast, we found decreased site-speci city in later time periods, especially in middle latencies (between 77 and 272 ms in the younger group, and between 87 and 285 ms in the older cohort). The cross-site similarity of these mid-latency potentials decreased signi cantly after removal of the AEP component, but remained mid-range and higher than both the earlier and later elements of the TEP.
By analyzing the residual mid-latency components, we found that there is prominent modulation of alpha band oscillation during this time period, which is present in both active and sham conditions. To our knowledge, alpha modulation exactly like this has not been described previously following non-M1 targets and sham stimulation, although there is some precedent for alpha and beta modulation with single pulse TMS to M1 [13-14,52, for reviews 18,71-73], which has been attributed to cortico-cortical neuronal processing [74][75][76], or to sensory return from muscle twitches [52,[77][78]. Because the midlatency components in both cohorts appear to synchronize after the TMS pulse but are not time-locked to the pulse, we describe this activity as non-speci c modulated alpha. Notably, this alpha modulation is similar across stimulation sites but distinct across subjects, suggesting that it re ects subject-speci c intrinsic oscillatory dynamics.
Alpha band modulation in the post-TMS period could be caused by some aspect of the stimulation that may not be transcranially evoked. The TEP recording in this window includes transcranial activity, sensory potentials, and attentional and behavioral responses to stimulation events. It is outside the scope of this paper to de ne all modulated oscillatory dynamics in the post-TMS period. However, brie y here we speculate about some explanations for the mid-latency components that were revealed with the presented AEP-removal method. One possibility is endogenous alpha phase resetting, due to subject speci c inhibitory responses to activation of cortical circuits, both transcranially and sensory-evoked. Suppression of cortical networks following TMS, re ected in the N100 time window, may be mediated through GABA-B receptors [79]. Other possible explanations include non-speci c modulation of arousal state, micro-blinks or micro-saccades [80][81][82], secondary downstream effects of auditory stimulation, predictive auditory processing of repetitive sounds, and other forms of event processing [83-86].

Theoretical Concerns
All in uences of auditory stimulation during an experiment cannot be isolated using these methods and need to be considered carefully in the context of experimental parameters and goals. For example, auditory stimulation can impact aspects of attention [86-87] and motor system excitability [88][89].
Further, if the design captures behavioral response, such as reaction time to a visual stimulus, intersensory facilitation should be accounted for [90][91][92][93]. At a neurophysiological level, the ICA technique can be used to suppress the AEP in the primary auditory region. However, the secondary downstream effects may persist (and indeed, may account for the observed alpha synchronization seen in middle latencies even in the sham). In the present study, we describe a method that can be used to examine in detail three dominant contributors to the TEP waveform (transcranial evoked potentials, sensory-evoked potentials, and site-independent modulation of ongoing oscillations), and not all consequences of auditory stimulation. The described technique can be used to study the three dominant parts of the TEP waveform with the caveat that this compartmentalized view of brain network activity may not be ecologically valid due to the multisensory nature of TMS.

Limitations
One limitation of our approach is that the AEP components were identi ed and labeled as AEP in a nonautomated way. This is a pervasive concern with using ICA for EEG data cleaning-the standard is to classify components by eye using component properties such as time course, topography and spectral pro le. This practice can be time consuming and susceptible to human bias. Many algorithms are useful in identifying speci c non-brain artifacts such as eye and muscle, although routinely mislabel other artifacts such as heartbeat. Because AEP is a non-artifact brain component, the component properties can be similar to other neural components and more di cult to distinguish from TMS-evoked neural components. However, unlike other components in the TEP, AEP is very stereotyped in time course and topography and detection can be further aided by the inclusion of a sham stimulation block.
With use of the above-described time course, topographical and condition in-speci c features, robust and stereotyped AEPs can easily be identi ed, as demonstrated in Figure 1. In order to expedite data processing and limit human rater bias, we suggest an automated AEP component identi cation is needed that classi es the AEP based on these stereotyped features [60].
A second potential limitation is that there could be a cost of analyzing the TEP with the AEP removed if condition, population, or individual participant relevant information is being carried in the auditory response. The auditory response may be relevant for questions related to auditory cortical excitability, plasticity and excitability gain control. However, because the components removed are likely to be predominantly AEP-related, they can be compared post-hoc across conditions or between populations or individual participants. Implications TMS is increasingly recognized as having signi cant potential utility in the neurophysiologic characterization of neurological disorders, as well as in the characterization of typical mechanisms of network connectivity and plasticity. However, due to the presence of sensory-evoked potentials such as the AEP, adequate experimental designs as well as appropriate and effective preprocessing of the TEP are critical. Consistent with other recent studies, our results suggest that early-latency TEPs are stimulation site speci c and largely unaffected by the AEP, and therefore could be optimal for characterization of evoked responses across groups or stimulation conditions. However, later latencies in the TEP are more at risk for AEP contamination, and therefore effective AEP suppression or isolation is needed. With adequate preprocessing, later components of the TEP show features that are stable across time, speci c to the individual, and may be relevant for cognitive mechanisms; future investigations are needed to further characterize the properties of these later components of the TEP. In addition, the factors that affect AEP topography, time-course and amplitude could be studied to better understand the condition, population, or individual participant information that is carried in the auditory response. Figure 2 Group level impact of removing AEP in the younger cohort (N=10). (A) DLPFC condition TEPs pre (black) and post (blue) removal of AEP components. Scalp topography pre (above) and post (below) at indicated latencies. (B) IPL condition TEPs/topographies pre (black/above) and post (purple/below) removal of AEP components. (C) Sham stimulation condition TEPs/topographies pre (black/above) and post (green/below) removal of AEP components. (D) Global mean eld potential (GMFA) area under the curve plots pre (darker) and post (lighter) AEP removal in three time-windows, de ned by sham GMFA peak activity pre-AEP removal.

Figure 3
Speci city of TEPs (cosine similarity) before and after removing AEPs in the younger cohort. (A) Betweensubject similarity index grids pre-and post-removal of AEP components. (B) Between-subject similarity index average (calculated from half of the grid, excluding the diagonal) pre-and post-removal (left), and within-subject stimulation site similarity (calculated as the average of the diagonal of IPL/DLPFC, sham/DLPFC, and sham/IPL grids; right). (C) Time windows used for within-subject site similarity in three windows (early: 14-77.24 ms, mid: 77.24-271.72 ms, late: 271.72-400 ms), de ned using the post-removal GMFA from DLPFC (darker) and IPL (lighter). (D) Within-subject stimulation site similarity pre and post removal in early, mid and late latency windows.   TEP/topography of all mid latency components contributing more than 28% variance in mid latency window in sham stimulation condition (N=16) from full group analysis.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.