{'fraction_non_alphanumeric': 0.03607241863171224, 'fraction_numerical': 0.022177350354221567, 'mean_word_length': 5.102015244857471, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Probing regional cortical excitability via input–output properties using transcranial magnetic stimulation and electroencephalography coupling

Abstract The modular organization of the cortex refers to subsets of highly interconnected nodes, sharing specific cytoarchitectural and dynamical properties. These properties condition the level of excitability of local pools of neurons. In this study, we described TMS evoked potentials (TEP) input–output properties to provide new insights into regional cortical excitability. We combined robotized TMS with EEG to disentangle region‐specific TEP from threshold to saturation and describe their oscillatory contents. Twenty‐two young healthy participants received robotized TMS pulses over the right primary motor cortex (M1), the right dorsolateral prefrontal cortex (DLPFC) and the right superior occipital lobe (SOL) at five stimulation intensities (40, 60, 80, 100, and 120% resting motor threshold) and one short‐interval intracortical inhibition condition during EEG recordings. Ten additional subjects underwent the same experiment with a realistic sham TMS procedure. The results revealed interregional differences in the TEPs input–output functions as well as in the responses to paired‐pulse conditioning protocols, when considering early local components (<80 ms). Each intensity in the three regions was associated with complex patterns of oscillatory activities. The quality of the regression of TEPs over stimulation intensity was used to derive a new readout for cortical excitability and dynamical properties, revealing lower excitability in the DLPFC, followed by SOL and M1. The realistic sham experiment confirmed that these early local components were not contaminated by multisensory stimulations. This study provides an entirely new analytic framework to characterize input–output relations throughout the cortex, paving the way to a more accurate definition of local cortical excitability.


| INTRODUCTION
The cerebral cortex presents a modular network organization allowing better robustness, adaptivity, and evolvability of network functions (Liao, Cao, Xia, & He, 2017;Meunier, Lambiotte, & Bullmore, 2010). Some characteristics of these brain modules are strongly modulated by topological and cytoarchitectural features, including the density of pyramidal cells (Fernández-Ruiz et al., 2013), the distribution of coactivated synapses, or the architectonic configuration of the cell populations (Kajikawa & Schroeder, 2011;Murakami & Okada, 2006). The same modular system exists in the vertical direction, with a columnar organization subdivided into different layers. Each layer contains a specific distribution of neuronal cell types and connections with other cortical and subcortical regions. The differences in lamination shape the input and output connectivity of neuronal populations and delineate distinct but interconnected functional cortical areas (He et al., 2009;Meunier, Lambiotte, Fornito, Ersche, & Bullmore, 2009). This modularity of functional brain networks suggests that, at the system level, discrete cortical regions or networks are associated with specific dynamical properties, as can be defined by their inputoutput properties. It is indeed likely that functionally relevant nodes of a network share common input-output properties, reflecting the aggregated architecture of the subsystems components. This might form the neural bases supporting the emergence of adaptive behaviors, including sensory, motor, and cognitive functions.
Here, we refer to input-output properties as the spectrum of modulations of a cortical area's activity to varying input levels, that could either be endogenous (from another cortical area or subcortical structure) or exogenous (using external stimulation, e.g. using TMS). In all biological systems, one would expect that input-output properties follow a few regular patterns, associated with a specific physiological or behavioral phenotype. For instance, some systems are broadly sensitive, that is, the output slightly changes, either linearly or nonlinearly, over a wide range of inputs, while other systems are ultrasensitive or bistable, that is, the output characteristics vary rapidly across a narrow range of inputs (for few examples within different size scales and inputs, see, e.g., (Desmurget & Sirigu, 2012;Trebaul et al., 2018) using direct cortical stimulation in humans; (Hegerl & Juckel, 1993) using peripheral stimulation in humans; (Yi, Wang, Tsang, Wei, & Deng, 2015) using single neuron model). This dynamical and networkspecific perspective on cortical physiology has been poorly or indirectly explored so far, but could provide a better estimate of the regional cortical excitability properties, characterizing the full inputoutput excitability profiles, from threshold to saturation.
The description of input-output properties of corticomotor or nonmotor neuron populations using TMS-EEG have been sparsely done though, and led to conflicting results. In an early study conducted by Komssi, Kähkönen, and Ilmoniemi (2004), the authors reported a nonlinear intensity dependency of the peak amplitudes of the overall brain response when stimulating the left and right motor cortices (but see (Saari et al., 2018) for opposite findings). Over the prefrontal cortex, Kähkönen, Komssi, Wilenius, and Ilmoniemi (2005b) found a linear dependency of the overall response on stimulus intensity, but with different peak latencies. In the two TMS-EEG studies, the authors reported similar potential distributions for the different intensities. In the time-frequency domain, it has been shown over the primary motor cortex (M1) that increasing TMS intensities induce a progressive synchronization of alpha and beta rhythm in both hemispheres (Fuggetta, Fiaschi, & Manganotti, 2005). Further data also showed that depending on the motor output (presence or absence of MEPs), the prestimulation EEG spectral power (Ferreri, Vecchio, Ponzo, Pasqualetti, & Rossini, 2014) and the interregional connectivity (Petrichella, Johnson, & He, 2017) differ. Additionally, this set of earlier papers revealed that evoked responses can be elicited even at subthreshold intensities (e.g., a minimal threshold of 60% resting motor threshold [rMT] was found by Komssi et al. to evoke a measurable brain activity over M1) although with different waveforms of overall activity (Kähkönen, Komssi, Wilenius, & Ilmoniemi, 2005a;Komssi et al., 2004;Komssi, Savolainen, Heiskala, & Kähkönen, 2007).
Paired-pulse TMS-EEG can also provide information on the local intracortical circuitry mediating inhibitory activity (Daskalakis et al., 2008;Farzan et al., 2009;Ferreri et al., 2011;Opie, Rogasch, Goldsworthy, Ridding, & Semmler, 2017;Rogasch, Daskalakis, & Fitzgerald, 2013;Rogasch, Daskalakis, & Fitzgerald, 2015;. These studies used a paired-pulse paradigm called long-interval intracortical inhibition (LICI) with two pulses separated by 100-200 ms to investigate the presumed activation of cortical GABAergic interneurons. For both the primary motor and prefrontal cortices, the mean cortical evoked activity was decreased and all typical components were found to be reduced compared to a single pulse TMS-EEG. By opposition, contrasting results have been published regarding the effects of another paired-pulse protocol called short-interval intracortical inhibition (SICI) over M1 using a 3 ms interval between the two pulses (Cash et al., 2017;Ferreri et al., 2011;Paus, Sipila, & Strafella, 2001). Although mediated by partially distinct receptors, a recent study showed similar amplitude reduction of the late TEP components induced by SICI and LICI (Premoli et al., 2018), which preclude any conclusion about the exact mechanisms of SICI over M1.
Here, we probed the activation of GABA A -ergic circuits over M1 by comparing it to single pulse TMS of different intensities and compared SICI modulations in three different brain regions.
Finally, a comprehensive definition of cortical excitability must consider the large interindividual variability in evoked neural responses, as increasingly reported in the literature (Gaspar, Rousselet, & Pernet, 2011). This variability undeniably limits the strength of the conclusions drawn from grand average ERP components, and motivates the use of more complex analytic tools (Bridwell et al., 2018), as we implemented in this article.
Hence, to fill the gap in the definition of cortical excitability, we provide a new dynamical and network perspective by characterizing the dynamical modes of the local source activity (LSA) evoked by TMS of increasing intensities in three distinct brain regions. We expected distinctive regional dynamical signatures with increasing intensities, potentially reflecting the recruitment of distinct neuronal populations. To consider the intersubject variability in ERP components, we also compared the quality of the linear regressions of the local evoked potentials on single trials. Its modulation through increasing stimulation intensities allows us to explore the sensitivity of the evoked neural activity across stimulation intensities and extract a new excitability threshold.

| Participants
Thirty healthy volunteers (19 males, aged 26.3 ± 6.2, two left handed) participated in the study. A first group of 22 subjects were recruited for the actual TMS-EEG experiment, while a second group of 10 subjects (including two participants from the first group) underwent a control experiment, in which a "realistic" sham procedure was used (see Conde et al., 2019;Gordon et al., 2018). All of them gave their written consent and filled a pre inclusion questionnaire screening for any contraindication for MRI nor TMS (Rossini et al., 2015). None had history of neurologic or psychiatric disorders, neither history of alcohol or substance abuse. All were free of any medicinal treatment likely to modulate their excitability. All participants received payment for their participation in the study. This study was approved by the ethical committee of Grenoble University Hospital (ID RCB: 2013-A01734-41), and registered on ClinicalTrials.gov (number NCT02168413).
The TMS-EEG experiment was performed in a 2-hr session. First, we prepared the subject for EEG (EEG cap setup) and we performed the coregistration step between the MRI and the physical subject's space that is necessary for the neuronavigation system. Second, a robotized hotspot hunting procedure was performed on a 7 × 7 grid (spaced by 7 m). The hotspot was defined as the cortical target maximizing muscular contractions from the FDI. The rMT was then assessed on the hotspot. Finally, we stimulated the three cortical areas at five different stimulation intensities (40, 60, 80, 100, 120% rMT, corresponding to 24 ± 4, 36 ± 6, 48 ± 8, 59 ± 10, 71 ± 11% of maximal stimulation output [MSO], respectively) and one paired pulse SICI protocol (conditioning pulse 80% rMT, stimulus pulse 120% rMT).
The order of stimulation sites was randomized across subjects and within each site, the order of the intensities was randomized but kept constant across sites. Participants had to sit still and relax during the measurements with their eyes open while staring at a black cross in front of them.
In order to take into account the last recommendations in the TMS-EEG field, an additional control experiment was run on 10 subjects (Belardinelli et al., 2019;Conde et al., 2019). None of them were naïve to TMS and TMS-EEG experiments: two of them underwent the actual TMS-EEG experiment described above, while the remaining eight took part in other on-going studies running in the lab.
The purpose of this control experiment was to provide a "realistic" sham, combining both auditory and somatosensory confounds (realSHAM). We electrically stimulated the right frontal area above right DLPFC, while placing the placebo coil on this latter target. Five different stimulation intensities were used and delivered using the same procedure previously described. rMTs were taken from previous TMS-EEG experiments undergone by these participants (40,60,80,100,120% rMT corresponding to 22 ± 5,33 ± 8,44 ± 9,55 ± 11,65 ± 14% MSO,respectively). An anecdotal evidence was found toward an absence of difference between the rMTs of the two groups (Bayesian independent t test, BF 10 = 0.6, see Section 2.11).

| TMS parameters
Biphasic TMS pulses induced an anteroposterior followed by posteroanterior current in the brain (AP-PA) using a MagPro Cool B65-RO butterfly coil (MagVenture A/S, Denmark) plugged in a MagPro ×100 TMS stimulator (MagVenture A/S). The coil was positioned and hold by a TMS robot (Axilum Robotics, France), navigated using Localite neuronavigation software (Localite GmbH). EMG electrodes were placed in a tendon-belly montage over the FDI. MEPs were recorded using a Dantec Keypoint portable EMG recording system (Natus Medical Inc.) and a CED micro 1401 MKII recording system (Digitimer, Cambridge Electronic Design, Cambridge, UK) for the six last subjects.
The rMT assessment was performed using the threshold hunting method (Awiszus, 2003) over the hotspot. The coil was positioned tangentially to the scalp surface in a posterior to anterior direction angled perpendicular to the central sulcus for M1. For the two other cortical targets, the coil was positioned perpendicular to the gyrus for DLPFC, and perpendicular to the axial plane for SOL. These angles were adjusted to the standard coil orientations used in the literature and to the mechanical constrains introduced by robot's motion (Janssen, Oostendorp, & Stegeman, 2015).
Each cortical target was stimulated at instantaneous frequency around 0.5-0.7 Hz for 2 min 30 s, resulting in an average number of 80-90 trials per stimulation point. For each cortical target, we adjusted the stimulation intensities using the Stokes formula calculating the scalp-cortex distance measured from subject's anatomical MRI (Stokes et al., 2007). We performed for each subject a classical sham condition (SHAM), which consisted of stimulations 3-5 cm over one of the cortical targets at the highest intensity used for each subject, allowing us to mimic the loudest sound generated by the TMS pulse without generating any somatosensory costimulation. We used active noise cancelation intraauricular earphones (Bose QC 20) combined with white noise to mask the TMS click susceptible to evoked auditory responses on the ongoing EEG activity [ter Braack et al., 2015]. The sound level was adjusted for each subject, so that the TMS click delivered at the loudest intensity (SICI) during the session became barely audible while the delivered sound was not loud enough to induce any discomfort. A thin layer of soft plastic was F I G U R E 1 Local TMS evoked potentials (TEPs) (a) and local source activity (LSA) power (b) for each stimulation site (lines) modulated by increased stimulation intensities (colors): from 40 to 120% resting motor threshold (rMT) and for 120% rMT and short-interval intracortical inhibition (SICI) on the left and right columns, respectively. The SICI condition corresponds to the response to the conditioned and test pulses (without substraction of/normalization to the test pulse only). Lines and shaded areas represent mean and SEM of local TEPs (a) and LSA power (b) z-scored against baseline. Black bars indicate periods of significant difference between conditions (see text). Bottom: Electrodes clusters taken for each site for the computing of local TEPs (left) and localization of the three scouts (regions of interest [ROIs]) defined for extracting LSA (right) in one representative participant placed on the coil surface to dampen both sensory and auditory feedbacks to the subject.

| Realistic sham stimulation parameters
The realistic sham stimulation was delivered using the MagPro Cool B65-A/P RO butterfly coil (MagVenture A/S), which is a coil originally designed to perform double-blind studies (Figure 3a). The coil was flipped on the placebo side. Concurrently to each TMS pulse, an electrical stimulation was delivered through two skin electrodes (stimulating area of 10 × 6 mm 2 ) placed on the scalp underneath the EEG cap above the DLPFC area, in a bipolar montage near electrodes AF4 and F6 (Figure 3a,b). Using this system, the electrical stimulation consists in a dissymmetric triangular monophasic pulse, with rise and fall times of 200 and 2,000 μs, respectively. The current intensity can be set from 0 to 6 mA, using a maximum voltage of 1,000 V. This intensity is adjustable by users (on an arbitrary scale from 0 to 10), and varies linearly together with the % MSO used in each experimental condition (40, 60, 80, 100, 120% rMT). The maximal intensity was defined for each subject prior to the EEG recording session. To that end, we first set the stimulator to 120% rMT, and then gradually increased the electrical stimulation intensity from 0 to 10. Subjects were asked to tell which stimulation intensity produced muscular twitches or skin sensations comparable in terms of strength, pain, or discomfort, to active TMS pulses. The selected intensity was then used throughout all the control experiment (mean 5.3 ± 3.9).

| EEG acquisition
EEG was recorded using a 64 channels TMS compatible system (BrainAmp DC amplifiers and BrainCap EEG cap, Brain Products GmbH, Germany). The EEG cap set up was done following the 10-20 standard system. Electrode impedances were adjusted and kept under 5 kΩ using conduction gel. The impedance levels were checked throughout the experiment and corrected if needed during breaks between conditions. The signal was recorded using DC mode, filtered at 500 Hz anti-aliasing low-pass filter and digitalized at 5 kHz sampling frequency. During the experiment, the Fz and Afz electrodes were used as reference and ground, respectively. Channel coordinates were individually assessed using the neuronavigation software at the end of the experiment.

| EEG preprocessing
EEG signals were processed using Fieldtrip (Oostenveld, Fries, Maris, & Schoffelen, 2011) and Brainstorm3 (Tadel, Baillet, Mosher, Pantazis, & Leahy, 2011) software, and other custom scripts written in MATLAB (The MathWorks Inc.). EEG signals were preprocessed semiautomatically based on the methodology described in Rogasch et al. (2014), for each condition (three targets, six conditions, and one sham) and each subject. First, the channels showing electrical noise (flat signal or peak-to-peak amplitude superior to 100 μV) spanning more than 15% of the trials were discarded from the analysis (on average, 1.4 ± 3.2 channels per condition). EEG signals were then epoched around the TMS pulse, using a −1 to +1 s time window of interest. TMS artifacts were discarded by cutting out the −5 to +17 ms period surrounding the TMS pulses. Two rounds of independent component analysis (ICA) were then applied in order to remove noise remaining in the signal. The first ICA suppressed the muscle artifacts, while the second ICA aimed at removing the decay artifact, ocular activity, auditoryevoked potentials, and other noise-related artifacts (Rogasch et al., 2014). Before the second ICA, the signal was spline interpolated over the −5 to +17 ms period, band-pass filtered (1-80 Hz), re-referenced using the average reference, and cleaned from bad trials (leading to a mean of 73.9 ± 9.7 trials left per condition). The ocular components were automatically identified using a threshold of 0.7 on the correlation product ρ between the spatial topographies of the components and a template of typical horizontal eye movements and blinks build from our own database by averaging over subjects. Other artifact components (decay, auditory-evoked potentials, and other noises) were detected by thresholding the z-score (above 4) of their mean activity against the prestimulus period, and by visual inspection. On average, 9 (±4.2) components were removed from the signal. Cleaned EEG time series were reconstructed using the remaining components and any isolated channel still showing remaining noise was discarded from further analysis. Time series of rejected channels were finally inferred using the activity averaged over their neighboring channels (see Harquel et al., 2016, figure 3, for an illustration of the main preprocessing steps).
Additionally, for the realistic sham data that are affected by a pronounced decay artifact, we applied a decay subtraction procedure between the two rounds of ICA (adapted from Conde et al., 2019).
Briefly, this procedure consists in subtracting the best fit of a twoexponential function from each trial of each channel. We used the nlinfit() function from MATLAB to estimate the five coefficients of the following regression function: with x being the time series of a specific trial and channel. Since the timing of the decay varies across conditions and channels, the fitting was optimized by processing it on increasing time window widths, from 200 to 800 ms by step of 100 ms. The width minimizing the mean squared error between the actual signal and the fitting function during the whole period of interest (0-1,000 ms) was taken.

| Global mean field potentials, TEPs, and LSA
First, to assess the TMS-evoked global cortical response, the global mean field potentials (GMFPs) were computed using the following formula: where t is time, C is the number of channels, V i is the voltage in channel i averaged across participants, and V mean is the mean of the voltage in all the channels.
Next, TEPs were computed for each target, stimulation intensity and subject by averaging the EEG signal across trials, using a baseline normalization (z-scoring) over the −200 to −5 ms period.
Grand average TEP was obtained by averaging normalized TEPs across subjects.
Source reconstruction for each nonnormalized TEP was performed following the default procedure proposed in Brainstorm 3 software (Tadel et al., 2011). First, the cortex and head meshes (15,000 and 10,000 vertices, respectively) of each individual were generated using the automated MRI segmentation routine of FreeSurfer (Reuter, Schmansky, Rosas, & Fischl, 2012). In order to extract LSA power, regions of interest (ROIs) were created on each individual anatomy using a mean spatial extent of 10 cm 2 , covering about 50-60 vertices of cortical mesh. LSA power was then computed for each cortical target by averaging the absolute, smoothed (using a spatial smoothing filter with full width at half maximum of 5 mm) and normalized (z-score against baseline) source activity within its corresponding ROI. Grand average LSA power was finally calculated for each stimulation site and intensity by averaging LSA power across subjects.

| LSA mode analysis
In order to disentangle the EEG response characteristics of various stimulation intensities through the identification of modes, we proceeded to a group ICA analysis over subjects for each stimulation site independently. Following the same methodology used in Harquel et al. (2016), this decomposition was performed on the signed LSA time series S i k of each stimulation intensity i and subject k, from −50 to +400 ms. The signed LSA time series were computed by averaging the signed and normalized source activity within each ROI (sign of sources with opposite directions were flipped before the averaging).
Each group ICA was performed after the concatenation of LSA matrices along the temporal dimension (Calhoun, Liu, & Adali, 2009) where N i is the number of intensities (6), N is the number of time bins (451), and K is the number of subjects (22). The matrix M was thus decomposed into N i (6) independent components (data dimension) using the logistic infomax ICA algorithm (Bell & Sejnowski, 1995) with the natural gradient feature from Amari, Cichocki, and Yang as implemented in EEGLab (Makeig, Bell, Jung, & Sejnowski, 1996).
Finally, the dynamical signature of each component was assessed in each individual by means of its time/frequency (TF) decomposition obtained using Morlet wavelet transform between 7 and 45 Hz (window width of 7 cycles, 0.5 Hz bandwidth). Individual TF power maps were normalized (z-score against baseline) and averaged across subjects.
2.9 | Linear regressions of early components of the local TEP in single trials Different linear regression analyses were performed at the scalp level.
First, the local TEPs x i were derived for each stimulation intensity i and each subject from the corresponding TEPs by averaging the signal of the four closest electrodes to each stimulation site ( Figure 1a).
For sham condition, local TEPs were extracted on central electrodes C1, Cz, and C2. The local TEPs were computed from +17 to +80 ms, in order to exclusively encompass the early components of the evoked activity. Then, linear regressions of the local TEPs were performed for each site on single trials s j extracted from the same electrodes and time window, so that: with (i, j) ∈ {40,60,80,100,120}% rMT, and then in a second analysis with (i, j) ∈ {80%,120%, SICI}, for each (i,j) intensity pairs. The term "paired intensities," used throughout this manuscript, refers to pairs where the intensity chosen to select a TEP matched the one used to select trials (see Figure 4c for a graphical description). In such cases, the TEP was thus computed from these same trials. For sham condition, the local TEP was regressed in its corresponding trials. Finally, the quality of the linear regression was assessed by extracting tstatistics associated with the local TEP x i factor, for each trial, intensity pair, site, and subject. For group analysis, these scores were averaged across trials for each intensity pair, site, and subject.

| Linearly scaled TEPs (simulated data)
We generated a set of simulated data whose components are linearly scaled with stimulation intensities in order to rule out a simple scaling effect of evoked temporal or spectral components. The set of simulated data ŝ i (t) were generated on a −400 to +600 ms period for each intensity (Figure 4b). The effect of the stimulation intensity consisted in a simple scaling of its inherent components' amplitudes, mimicking what is usually observed when increasing intensities on sensory evoked potentials (Juckel, Csépe, Molnár, Hegerl, & Karmos, 1996;Shiga et al., 2016;Tsuji, Lüders, Dinner, Lesser, & Klem, 1984). Waveform of simulated signals was designed to get close of what is typically reported in TEPs (Farzan et al., 2016), that is, from two to six alternative components together with some oscillatory patterns. Since the sole aim of these simulated data was to test the amplitude-scaling hypothesis mentioned above, the level of complexity of our simulation was kept rather low. Simulated signals were composed of two evoked early components (positive and negative peaks p 1 and p 2 ) within the first 80 ms and a late induced oscillatory activity o, contaminated with noise ε drawn for a uniform distribution filtered in the 1-80 Hz frequency band: The two early components were modeled using Gaussian functions of different mean (30 and 45 ms) and SD (5 and 20 ms, respectively) parameters. The late induced activity consisted in a sinusoidal function mimicking an alpha rhythm oscillation (10 Hz) starting from +150 to +350 ms. Finally, the effect of the stimulation intensity was modeled using the amplitude factor A i . Five different intensities were modeled, where A i ∈ {20,40,60,80,100} (Figure 4b).

| Statistics
Statistical analyses were conducted using the Fieldtrip and MATLAB statistical toolboxes on EEG signal, and using JASP Team (2018) (Version 0.9) for the Bayesian statistics analysis of regression quality scores.

| Local source activity
For each stimulation site, significant differences in the LSA across stimulation intensities were assessed over time, from +17 to +400 ms, using nonparametric permutation tests. The effect at the sample level was evaluated using the dependent samples F-statistics and T-statistics for the comparison of all stimulation intensities and the comparison between the SICI and 120% condition, respectively. The significance probability was then inferred using Monte-Carlo procedure with 10,000 permutations.
Finally, p-values were temporally corrected for multiple comparisons: differences were considered as significant at p < .05 for at least 20 consecutive time bins (20 ms, see Blair & Karniski, 1993;Carota et al., 2010;Harquel et al., 2016). Statistical significance of TF maps of ICA components was obtained using paired comparisons against baseline. A nonparametric Wilcoxon test was performed per time-frequency bin, and the resulting p-values were spatiotemporally corrected: differences were considered significant for p < .05 for at least five consecutive frequency bins and 20 time bins (tiles of 2.5 Hz × 20 ms).

| Regression quality scores
The mean regression quality scores across trials were analyzed using the Bayesian equivalent of repeated measures analysis of variance (rmANOVA) and analysis of variance (ANOVA) tests. Additional post hoc analysis was performed using the Bayesian equivalent of independent and paired t tests. Three analyses were conducted. The first one took all data expect SHAM and SICI conditions as inputs, and performed an ANOVA with three fixed factors: TEPs' stimulation intensity, single trials' stimulation intensity, and stimulation site (including active sites and realSHAM). Subjects were included as a random factor. Then, a second ANOVA analysis focused on the regression scores obtained in paired intensities (e.g., when the TEP of 80% rMT was regressed in its corresponding trials of 80% rMT, see above) using all data (including SHAM or SICI condition), with two fixed factors: stimulation intensity and stimulation site (including active sites and realSHAM). Subjects were included as a random factor as well. Finally, a third analysis was conducted specifically on SICI, 80 and 120% rMT conditions, using a rmANOVA with three factors: TEPs' stimulation intensity, single trials' stimulation intensity and stim-

| RESULTS
Overall, the 22 participants tolerated well the experiment. For two subjects, however, the SICI condition over the DLPFC was too painful, and was therefore omitted. The 10 additional subjects recruited for the realistic sham experiment did not report any adverse effect. However, three of them reported an increase of the pain throughout the procedure, probably due to the cumulative effect of the electrical stimulation on the skin. Below, different aspects of the input/output properties are described through complementary features computed for each targeted cortical area.  The maximal voltage amplitude of early components (17-80 ms, across all electrodes) in active conditions were 2.2 ± 0.7, 2.5 ± 0.9, 3.4 ± 1.7, 4.3 ± 2.1, 6.4 ± 5.3, 7.8 ± 6.5, and 6.8 ± 3.8 μV for SHAM, 40%, 60%, 80%, 100%, 120%, and SICI condition, respectively. In contrast, the maximal voltage amplitude of early components (17 to 80 ms, across all electrodes) in realistic sham conditions were 1.4 ± 0.3, 1.8 ± 0.5, 1.8 ± 0.5, 2.1 ± 0.6, and 2.5 ± 0.8 μV for 40, 60, 80, 100, and 120% stimulation intensities. Finally, the distribution of the electrical fields differed across stimulation sites, both in terms of spatial and temporal features.

| LSA power
Sources of TEPs for the three stimulated regions were estimated and local cortical responses (LSA) were extracted from the mean source time series of an ROI, centered on the stimulation target (Harquel et al., 2016). Figure 1b shows  The LSA response patterns to SICI also differed in the three regions (right part of Figure 3b). In M1, late activities (100 and 300 ms) were significantly inhibited by SICI (p < .05, t test corrected for multiple comparisons) as compared with the LSA over the same region obtained with a stimulation intensity of 120% rMT. In contrast,  In contrast, we found that the variance of simulated data was mainly First, an extreme effect of site (BF incl > 10 13 ) was found, suggesting that the regression fit is different in the four regions ( Figure 6a). Post hoc comparisons showed with strong evidence that M1 and SOL had similar fits (BF 10 = 0.07), which both exceeded DLPFC (vs. M1: BF 10 > 500, vs. SOL: BF 10 > 10 3 ). All active conditions showed stronger regression quality scores than realistic sham (BF10 = 9.8; >10 5 ; >10 7 for DLPFC, M1, and SOL, respectively). The two other main factors (intensity of TEPs and trials) were also significant (BF incl > 10 13 for both). This confirmed that the quality of the regression differed with intensity for both the TEP and the trials used for the regression, as the signal to noise ratio (SNR) gradually increased with intensity in the EEG signal.

The significant TEPs intensity by trials intensity interaction
(BF incl > 10 13 ) indicated that the regression of a given TEP fits better with the trials corresponding to the same intensity (paired intensities). Moreover, for all three active sites together, post hoc comparisons showed that the quality of regression was maximal for paired intensities (see Figure 6a), except for 40% rMT. Above 40% rMT, the best regression qualities were systematically obtained when using the same trials intensity than the TEPs intensity used for the regression (with strong to extreme evidence), confirming that each intensity has its own spatiotemporal signature. This was not the case for realistic sham conditions, where this one-to-one association was only observed in the 60 and 120% conditions (with moderate to strong evidence). Finally, a triple interaction TEPs intensity by trials intensity by site interaction (BF incl > 10 4 ) showed with extreme evidence that the increase in stimulation intensity had a different impact on the regression fit in the four regions.
We conducted a complementary analysis on the regression quality scores obtained using the late components (from 80 to 400 ms) of the central evoked potentials, on C1 Cz and C2 electrodes (supplementary Figure S1). A similar Bayesian ANOVA was performed, that showed extreme evidence (BF incl > 10 13 ) for all the three main effects and the interaction between TEPs intensity and trials intensity. However, moderate to strong evidence suggested that the interaction between site and TEPs intensity (BF incl = 0.17), as well as the triple interaction (BF incl = 0.07), had no effect on the model. between TEPs and trials intensity did not reveal any systematic better fit on paired intensities. The best fit was obtained for trials intensities that were equal to or higher than the TEP used (see Figure S1, Supporting Information). For each TEP above 40% rMT, we found extreme evidence toward differences between paired and lower intensities (all BF 10 > 100), while a lack of evidence or a moderate evidence toward absences of difference emerged from the comparison between paired and higher intensities (BF 10 between 0.17 and 1).

| Regression quality scores on simulated data
To demonstrate that increased intensities do not act like a simple scaling of the evoked components on real data, we generated a set of simulated data whose components were linearly scaled with stimulation intensities (see Section 2). The same statistical model applied to the simulated data showed significant main effects of TEPs intensity (BF incl = +inf) and trials intensity (BF incl = +inf) as well as a significant TEPs intensity by trials intensity interaction (BF incl > 10 7 ). Importantly, post hoc tests revealed that the best fit was obtained with 120% rMT trials whichever TEPs intensity was used for the regression. Figure 6a shows this clear linear relationship between intensities and regression quality scores, and the difference regarding curve shape with real data.

| Regression quality scores in paired intensity
Next, we restricted our model to the TEPs regressed with their corresponding trials (e.g., M1 TEP 100% rMT regressed on M1 100% rMT trials), including the sham and the realistic sham conditions in the model. A Bayesian ANOVA revealed extreme main effects of TEPs intensity (BF incl > 10 14 ), site (BF incl > 10 11 ), and a TEPs intensity by site interaction (BF incl > 10 5 ). Figure 6b shows the four different response curves associated with the four stimulation sites. Post hoc comparisons showed that a stimulation intensity effect was present for all three active sites (for M1: BF 10 > 10 12 , for DLPFC: BF incl > 10 7 , for SOL: BF incl > 10 5 ), The analysis of Figure 6b suggests a saturation effect at 100% rMT for DLPFC and SOL. Post hoc pairwise comparisons indeed provided substantial evidence for similar regression fits between 100 and 120% over DLPFC (BF 10 = 0.23) and SOL (BF 10 = 0.23), and moderate evidence for a better regression at 120% for M1 (BF 10 = 6.5). They also revealed that the minimal intensity needed to reach a statistical difference with sham was 60% rMT over M1 (Figure 6c). Over the DLPFC and SOL, TMS needed to be applied at 100 and at 80% rMT, respectively, to reach strong statistical evidence. No statistical evidence was found to infer about a difference or an absence of difference between realistic sham conditions and sham, in any used intensity (all BF 10 fell between 0.38 and 0.74, Figure 6c).

| Regression quality scores in SICI condition
The last rmANOVA performed on 80, 120% rMT and SICI conditions also indicated a significant TEPs intensity effect (BF incl = +inf) and a significant site × TEPs intensity interaction (BF incl = 9.2). Post hoc tests (Figure 7a) showed that the regression quality scores acquired with SICI were smaller than for 120% rMT only for M1 (BF 10 = 10.6).
Moderate evidence tended to show that it had no effect on the quality of the regression for DLPFC and SOL (BF 10 = 0.25 and BF 10 = 0.23, respectively). Considering all intensity pairs, a last analysis showed again a typical response pattern for the SICI trials (extreme TEPs intensity by trials intensity interaction: BF incl > 10 15 ) that were better regressed by their own TEP and could not be found in the other conditions ( Figure 7b). This comparison further suggests that SICI, which is composed by a first TMS pulse at 80% rMT and a second one at 120% rMT, induced a specific pattern of activity which significantly differed from single pulse TMS given at 80 or 120% rMT.

| DISCUSSION
In the present paper, we reported a set of new EEG markers able to quantify interregional differences in input-output properties and inhibitory activities in three different brain areas.
F I G U R E 7 Regression quality scores modulated by short-interval intracortical inhibition (SICI). (a) Regression quality scores obtained for each site in paired intensities. (b) Regression quality scores obtained for each site, trials intensity (y axis), and TMS evoked potentials (TEPs) intensity (panel columns). In each panel column, the corresponding paired intensities (where i = j, see Section 2) are highlighted using its specific color 4.1 | Regional input output properties of the cerebral cortex

| Local neural activity scales with stimulation intensity
LSA recorded in the three regions with the five intensities showed an overall dose-dependency relationship. This finding replicates previous TMS-EEG or TMS-fMRI studies, either stimulating the primary motor cortex (Bohning et al., 1999;Komssi et al., 2004;Komssi et al., 2007) or the left middle frontal gyrus (Kähkönen et al., 2005a). The earlier TMS-EEG studies have intrinsic limitations in terms of data acquisition and analysis, precluding a full understanding of the input-output response function to TMS. Although ROI analyses in the source space have been criticized (Farahibozorg et al., 2018), the LSA profiles computed from each of the three ROIs encompassing the stimulated area, allowed us to distinguish complex shapes of biological input-output relationships. Unlike Komssi et al. (2004) where all EEG peaks depended nonlinearly on stimulation intensity up to 100 ms poststimulation, our patterns of changes yielded to mixed results in M1. The N100 showed a clear linear relationship while the P30 increased nonlinearly with increased stimulation intensity. This discrepancy was also present in the two other regions showing nonlinear dependencies (i.e., the N100/P200 for both the DLPFC and SOL), as well as linear dependencies (i.e., the P200 and P30 for the DLPFC and SOL, respectively).
While each EEG component reflects brain activation and is thought to be associated with one or more cognitive processes (Brandeis & Lehmann, 1986;Sur & Sinha, 2009) Comparing M1 and DLPFC, Kähkönen et al. (2005b) found different reactivities of motor and prefrontal cortices and different dynamics. A quadratic polynomial function described the data of motor cortex, whereas a linear model fitted better for the response-stimulus intensity function of prefrontal TMS. We also found that the dynamic of the input/output functions differed across the three regions, suggesting different dynamical properties of neuronal responses. This provides additional evidence that TMS-EEG can noninvasively probe regional differences in cortical microcircuits underlying functional cytoarchitecture (Harquel et al., 2016). Of interest, a recent study investigated the effect of stimulation intensities of intermittent theta burst stimulation on the cortical properties assessed with TMS-EEG.
The authors reported an inverse U-shaped relation between intensity and induced plastic effects, where 75% iTBS yielded the largest neurophysiological changes (Chung et al., 2018). These results not only raised interesting aspects about the relationship between intensity and plasticity induction but also about homeostatic regulation maintained through the recruitment of excitatory and inhibitory subpopulations of neurons, which can be characterized by TMS-EEG components (Premoli et al., 2014). Then, for rTMS treatments applied in "silent" regions, a systematic description of the input-output functions of different cortical areas is crucial in order to induce the most efficient plasticity change. An individual and fine-tuning of rTMS intensity is particularly important given that the brain responses to rTMS are dose dependent, precisely, stimulation intensities influences the plasticity induction (see, e.g. Nettekoven et al., 2014or Fitzgerald et al., 2002). Using regional biophysical models of neural plasticity induced by TMS and individual input-output excitability profiles, getting more reliable rTMS outcomes might be possible. Alternatively, for clinical practices, the use of region-specific atlases of excitability profiles might already help the clinicians to define the optimal ranges of stimulation intensities, although this last option does not take into account the particularities of pathological brains. that each brain region mostly resonates at its own natural frequency.

| Intensity-dependent spectral properties of the evoked neuronal responses
Moreover, our previous work showed that distant cortical areas can also share common dynamical properties, depending on their local cytoarchitectonics (Harquel et al., 2016). In the primary motor cortex, Fuggetta et al. (2005) showed that different stimulation intensities appear to involve different levels of modulation of oscillatory activity.
Precisely, increases in alpha and beta power were found to be more pronounced with increasing of stimulation intensity from subthreshold to the 130% rMT. Here, we used LSA modes decomposition to further reveal that different intensities in each of the tested regions were associated with distinct and complex patterns of oscillatory activities, rather than with a linear or nonlinear scaling of a specific frequency pattern.
For the three tested sites, our group ICA including the five intensities and the SICI condition revealed at least four components with comparable level of explained variances, arguing for the existence of mixed and complex dynamics. Moreover, each component was mainly driven by one unique stimulation intensity, which entails its own oscillatory signature. In addition, to explain approximately 90% of the total variance, all components have to be included, suggesting that each stimulation condition is associated with a complex mixture of oscillatory signals.
As a matter of fact, the power of the so-called "natural" rhythm of a given area was not the only one impacted by stimulation intensity.
Instead, this analysis showed that the input-output relationship is unlikely to be driven by a simple linear scaling of the stimulation intensity as for the simulated data. The variance of simulated data was mainly explained by one single component (explained variance on the simulated data: 90.9, 3.5, 2.7, 1.2,and 0.78%), which means that the same information content, and its related oscillatory signature, can be found in all conditions. To sum up, this analysis demonstrated that each stimulation intensity induced oscillatory activities reflecting complex combination of frequency bands, distinct from each other and different across brain regions.
The changes observed in the EEG activity at low TMS intensity are thought to originate from the stimulation of the superficial layers of the cortex through both direct and indirect excitation of pyramidal neurons in the gray matter (Di Lazzaro et al., 2000;Kujirai et al., 1993;Ziemann, Rothwell, & Ridding, 1996). In contrast, by increasing stimulation intensities, the direct axonal pathways in deep gray matter structures get activated (Amassian & Cracco, 1987;Nakamura, Kitagawa, Kawaguchi, & Tsuji, 1996). This could additionally activate subject. This method allows to fully consider the intersubject variability of the dynamics of the evoked response, since subject-specific dynamical contents of TEPs will not influence the regression process.
The quality of the regression could provide a new local readout for cortical excitability and dynamical properties specificity. Precisely, cortical excitability could be defined by the quality of the regression of TEPs' early components compared across intensities for each cortical region. The rationale of this new metric is that, at similar stimulation intensities, highly excitable neural populations would be more prone to produce electrical activity above noise level in a single trial basis, compared to low-excitable populations. This will be associated with better quality of TEP regression. In the same line, inspecting the shape of the relationship between regression quality and stimulation intensity might provide a more specific definition of cortical excitability definition.
The results showed first that regression fits were overall smaller in the DLPFC compared to SOL and M1, indicating a lower excitability in this region, at least regarding early components. This means that each trial carries more information from local neural activations for M1 and SOL compared to DLPFC in which a single trial poorly explains the average TEP. A new regional cortical excitability index can then be inferred based on the comparison of the regression quality obtained in the sham condition and with different TMS intensities.
This index would correspond to the weakest stimulation intensity able to elicit significant TEP regression quality (i.e., different from sham stimulation). Note that as expected, our analysis also showed that the goodness of fit gradually improves with intensity, reflecting a better SNR with increasing intensities. This metric has the advantage to be defined locally even in "silent" areas and does not rely on peripheral readout. It is biologically informative because it directly reflects the quantity of energy needed to evoke a meaningful EEG signal in a given brain region. This idea is reinforced by the fact that no significant information could have been retrieved from realistic sham conditions. As shown by the comparison between areas, it is sensitive enough to discriminate the different excitability levels between different cortical areas. The input-output curves drawn from regression score suggested that the primary motor cortex was the most excitable area (eliciting significant activity from 60% rMT), followed by SOL and DLPFC (80 and 100% rMT, respectively). Fecchio and colleagues also report larger local mean field potentials evoked over M1 compared to prefrontal, premotor, and parietal targets (Fecchio et al., 2017). This has direct consequences on TMS titration for clinical trials. Our data suggest that adjusting rTMS intensities (for instance, applied to the DLPFC) to the rMT is suboptimal. The new excitability metric we present here brings instead an accurate estimate of regional excitability, which could serve as a basis to better adjust stimulation intensities. In the same vein, Casali et al., 2010 also derived a local excitability estimate based on the minimal TMS intensity needed to significantly activate more than 1% of the cortical sources within the stimulated area, and found comparable values over the superior occipital lobe (Casali et al., 2010).
Furthermore, regression quality could also provide interesting insights about the dynamical properties of the evoked response and especially about its modulation across brain areas and stimulation intensities. Our results showed that above 40% rMT, the quality of regression is maximal for the trials corresponding to the regressed TEP, confirming that each intensity has its own spatiotemporal signa- 4.3 | Regional inhibitory properties of the cerebral cortex 4.3.1 | Local neural activity shows opposite response patterns to SICI Paired-pulse TMS delivered with an interval of 2-3 ms can noninvasively probe the level of GABAA receptor (GABAAR) mediated inhibition (Di Lazzaro et al., 2007;Kujirai et al., 1993). Importantly, the resulting EEG responses reflect the combination of the conditioning effects of the first pulse on the second one and vice versa. Then the interpretation is not as straight forward than a unique activation of GABAAR.
To our knowledge, there are only four studies investigating the EEG correlates of SICI in M1 (Cash et al., 2017;Ferreri et al., 2011;Paus et al., 2001;Premoli et al., 2018) and one in the DLPFC (Cash et al., 2017 (Premoli et al., 2018), similar to earlier findings related to LICI (Premoli et al., 2014). Because SICI has distinct impacts depending on the stimulated area, TMS-EEG could provide a new regional readout for drug testing specifically targeting GABAergic mechanisms in predefined brain areas.

| Potential experimental confounds
Our comparison of the EEG effects induced by the different intensities is inevitably confounded by TMS-induced muscular, auditory and somatosensory responses, as called peripheral evoked potentials (PEPs). Below, we will discuss how our data and our new analysis approach might add new insights concerning the contribution of these multisensory costimulations into the TMS-induced EEG response, with respect to the classical sham procedure and the realistic sham experiment we conducted. Several recent lines of evidence suggest that realistic sham stimulation induces a cortical response pattern close to the one evoked by real TMS over the scalp Conde et al., 2019]. The current debate on this topic (Belardinelli et al., 2019;Siebner et al., 2019) prompted us to better disentangle the multisensory temporal and spatial response patterns from the real transcranial evoked brain response to TMS.
First, the TMS-induced auditory sound may have contributed to the TMS intensity-dependent changes of components' amplitude in the three tested regions, contaminating the early components (i.e., P30, P60) and the late components (i.e., the N100 and P200).
Indeed, the "click" sound of TMS increases with intensity (Dhamne et al., 2014) producing an increase in AEP (Juckel et al., 1996). Additionally, intensity-dependent artifact can result from cranial muscle activity related to direct depolarization of muscle fibers by the TMS pulse or from activation of the nerves innervating the muscles (Mutanen, Mäki, & Ilmoniemi, 2013). Finally, there are additional regional sources of artifacts. For example, over M1, the presence of intensity-dependent muscle reafferent inputs to S1 induced by suprathreshold TMS over M1 results in sensory-evoked potentials in the EEG, which can contaminate the TEP (Fecchio et al., 2017). In the same line, TMS over the DLPFC can be uncomfortable and this feeling of discomfort is proportional to TMS intensity. One way to quantify the weight of these nonneuronal signals in the various TEPs would be to systematically rate the discomfort and pain induced by all the stimulation conditions and relate it to the associated TEPs.
A further development of our cortical mapping would be to systematically and online fine-tune TMS parameters (angle, intensity, position) to minimize artifacts and maximize cortical responses before starting the acquisition [see, e.g., Casarotto et al., 2016]. This approach would allow a better comparison of interregional signals and a less extensive use of postprocessing computations. However, we tried in the present work to minimize confound effects both during acquisition and EEG processing steps. First, we used a realistic sham condition consisting in a TMS sham coil able to mimic the TMS multisensory effects without directly stimulating the brain (Smith & Peterchev, 2018). We also used state-of-the-art methods to reduce the auditory component (plastic form under the coil and sound-protective headphones playing white noise). Then, we applied state-of-the-art methods to process TMS-EEG data by using the methodology developed in 2014 by Rogasch et al., which rely on systematic and meticulous data cleaning steps using two ICA rounds. Such a method allowed us to remove any residual muscular or auditory artifacts from our data, and to compute clean TEPs showing significant early neuronal activations (Belardinelli et al., 2019).
Despite all these measures, it is still impossible to dissociate the multiple sources of multisensory stimulations induced by TMS. Then the resulting input-output patterns are partly ambiguous and cannot be definitively attributed to direct local cortical TEP profiles. However, the results of our realistic sham experiment clearly stated that no significant information can be drawn locally from the early components of PEPs (<80 ms) using our methodology, contrary to its late central components.
The quality of regression of the realistic sham conditions on early components never significantly differed from noise, as quantified with the classical sham procedure used in the first experiment. This is in line with recent findings suggesting that only late components appear to contain significant PEPs Freedberg et al., 2020). In contrast, recruitment curves drawn from active stimulation differed from noise and showed different patterns across sites, possibly revealing different input-output properties of the cortical tissue. Additionally, both our LSA modes and regression quality analyses enabled us to link one stimulation intensity with one specific dynamical signature, distinct from each other, which was not the case with realistic sham data. Using simulated data, we also demonstrated that these results were in contradiction with the hypothesis of a simple scaling effect of evoked temporal or spectral components, which would be the case for increasing auditory or somatosensory stimulation (Juckel et al., 1996;Shiga et al., 2016;Tsuji et al., 1984). Interestingly, this latter hypothesis is partially confirmed on late central components, since our complementary analysis revealed that no specific dynamic signature can be drawn across (or between) stimulation intensities at such latencies. However, further analyses of these components have to be performed. Several pieces of evidence, such as the different spatial and temporal features of the late electrical fields across conditions, the main effect of site for the regression quality score, and the presence of late components induced by intracranial stimulations (e.g., Keller et al., 2018;Kunieda et al., 2015) support the idea that relevant information might be contained at such latencies.

| CONCLUSION
In this article, we examined the TMS intensity dependent effects on the EEG signals from a temporal and spectral perspective, and used a new analytic approach to derive regional input-output profiles of the EEG responses to TMS. We reported complex dynamical responses in three distant regions in terms of components' amplitude and oscillatory signatures. These complex properties of the local neuronal responses to TMS largely depend on the intrinsic cytoarchitecture and connectivity patterns of the stimulated area. Then, the systematic description of their intrinsic dynamical properties brings important knowledge into cortical physiology. Furthermore, our data have implication in clinical research. This new assessment of regional cortical excitability will help resolving new challenges especially in the context of pharmacological or brain stimulation induced modulations of cortical excitability (Devergnas & Wichmann, 2011;Karabanov et al., 2015;Nitsche et al., 2003). The primary motor cortex has been largely used as the experimental model to study brain reactivity to TMS and to normalize stimulation parameters when targeting nonmotor areas.
Our data show that the local EEG response may be highly specific, and should not be extrapolated to other brain regions. Then using the rMT as a basis to define stimulation intensities in rTMS protocols might contribute to the large variability reported by clinical rTMS trials. Based on linear regressions, a realistic regional threshold could be derived and could serve as a basis for to set stimulation parameters at the group level or even on an individual basis. Hence, the development of more accurate control strategies of TMS-induced changes in cortical excitability will eventually facilitate predicting the effect of rTMS applied to nonmotor brain areas.

ACKNOWLEDGMENTS
This work was funded by the Agence Nationale pour la Recherche grant "ANR-15-CE37-0015-01" and by NeuroCoG IDEX UGA in the framework of the "Investissements d'avenir" program (ANR-15-IDEX-02). Data were acquired on a platform of France Life Imaging Network partly funded by the grant "ANR-11-INBS-0006."