Intensity and Dose of Neuromuscular Electrical Stimulation Influence Sensorimotor Cortical Excitability

fnins-14-593360 January 15, 2021 Time: 14:6 # 1 ORIGINAL RESEARCH published: 15 January 2021 doi: 10.3389/fnins.2020.593360 Edited by: Julia Stephen, Mind Research Network (MRN), United States Reviewed by: Sidney Grosprêtre, University of Franche-Comté, France Dipanjan Roy, National Brain Research Centre (NBRC), India *Correspondence: Ainhoa Insausti-Delgado [email protected] Specialty section: This article was submitted to Brain Imaging Methods, a section of the journal Frontiers in Neuroscience Received: 10 August 2020 Accepted: 30 November 2020 Published: 15 January 2021 Citation: Insausti-Delgado A, López-Larraz E, Omedes J and Ramos-Murguialday A (2021) Intensity and Dose of Neuromuscular Electrical Stimulation Influence Sensorimotor Cortical Excitability. Front. Neurosci. 14:593360. doi: 10.3389/fnins.2020.593360 Intensity and Dose of Neuromuscular Electrical Stimulation Influence Sensorimotor Cortical Excitability Ainhoa Insausti-Delgado1,2,3*, Eduardo López-Larraz1,4, Jason Omedes5,6 and Ander Ramos-Murguialday1,7 1Institute of Medical Psychology and Behavioral Neurobiology, University of Tübingen, Tübingen, Germany, 2International Max Planck Research School (IMPRS) for Cognitive and Systems Neuroscience, Tübingen, Germany, 3IKERBASQUE, Basque Foundation for Science, Bilbao, Spain, 4Bitbrain, Zaragoza, Spain, 5Instituto de Investigación en Ingeniería de Aragón (I3A), Zaragoza, Spain, 6Departamento de Informática e Ingeniería de Sistemas (DIIS), University of Zaragoza, Zaragoza, Spain, 7Neurotechnology Laboratory, TECNALIA, Basque Research and Technology Alliance (BRTA), Donostia-San Sebastián, Spain Neuromuscular electrical stimulation (NMES) of the nervous system has been extensively used in neurorehabilitation due to its capacity to engage the muscle fibers, improving muscle tone, and the neural pathways, sending afferent volleys toward the brain. Although different neuroimaging tools suggested the capability of NMES to regulate the excitability of sensorimotor cortex and corticospinal circuits, how the intensity and dose of NMES can neuromodulate the brain oscillatory activity measured with electroencephalography (EEG) is still unknown to date. We quantified the effect of NMES parameters on brain oscillatory activity of 12 healthy participants who underwent stimulation of wrist extensors during rest. Three different NMES intensities were included, two below and one above the individual motor threshold, fixing the stimulation frequency to 35 Hz and the pulse width to 300 µs. Firstly, we efficiently removed stimulation artifacts from the EEG recordings. Secondly, we analyzed the effect of amplitude and dose on the sensorimotor oscillatory activity. On the one hand, we observed a significant NMES intensity-dependent modulation of brain activity, demonstrating the direct effect of afferent receptor recruitment. On the other hand, we described a significant NMES intensity-dependent dose-effect on sensorimotor activity modulation over time, with below-motor-threshold intensities causing cortical inhibition and above-motor-threshold intensities causing cortical facilitation. Our results highlight the relevance of intensity and dose of NMES, and show that these parameters can influence the recruitment of the sensorimotor pathways from the muscle to the brain, which should be carefully considered for the design of novel neuromodulation interventions based on NMES. Keywords: neuromuscular electrical stimulation (NMES), electroencephalography (EEG), afferent cortical activation, sensorimotor oscillatory rhythm, artifact removal Frontiers in Neuroscience | www.frontiersin.org 1January 2021 | Volume 14 | Article 593360 fnins-14-593360 January 15, 2021 Time: 14:6 # 2 Insausti-Delgado et al. NMES Parameters Influence Cortical Excitability INTRODUCTION Neuromuscular electrical stimulation (NMES) is an electrophysiological technique that consists of applying electrical currents on the skin to depolarize motor and sensory nerves beneath the stimulating electrodes (Bergquist et al., 2011). NMES has been used as a neuroscientific tool to study sensorimotor neural mechanisms from the muscles or peripheral sensory receptors (mechanoreceptors, nociceptors, etc.) to the spine and brain (Collins, 2007;Carson and Buick, 2020). NMES has also been utilized as a neurorehabilitative tool to reduce muscle atrophy, and to improve muscle tone and motor function in patients with paralysis after stroke (Knutson et al., 2015;Yang et al., 2019) or spinal cord injury (SCI) (Patil et al., 2015). The working principle of rehabilitative NMES is based on the following: (1) the direct effect on muscle tone and (2) the activation of receptors and sensory axons that send afferent volleys to the sensorimotor cortex, after being processed by spinal networks and subcortical structures. Numerous studies using functional magnetic resonance imaging (fMRI) and near-infrared spectroscopy (NIRS) have investigated the activity produced by NMES in different brain structures (Kampe et al., 2000;Blickenstorfer et al., 2009;IftimeNielsen et al., 2012;Schürholz et al., 2012;Wegrzyk et al., 2017). These works demonstrated that the brain activity is proportionally increased with the applied stimulation intensity (Backes et al., 2000;Smith et al., 2003;Schürholz et al., 2012). Recent experiments have also shown that peripheral stimulation can modulate corticospinal excitability, measured by motor evoked potentials (MEP) (Chipchase et al., 2011; Veldman et al., 2014). As the stimulation intensity increases, there is a progressive recruitment of more afferent receptors (i.e., cutaneous mechanoreceptors, muscle spindles, and Golgi tendon organs) that modulate spinal and cortical circuits to a different extent (Maffiuletti et al., 2008;Bergquist et al., 2011; Golaszewski et al., 2012). It has been proven that the afference provided by muscle spindles and Golgi tendon organs due to muscle contraction travels through the spinal cord to the somatosensory cortex and can directly project to the motor cortex (Carson and Buick, 2020). Therefore, the presence or absence of muscle contraction elicited by NMES has a direct impact on somatosensory cortex and, indirectly, on motor cortex excitability (Sasaki et al., 2017). The neural excitability depends on the modulation of nervous structures, such as spinal networks, involved in the afferent transmission from the stimulated muscle to the brain. Although the effect of intensity on cortical activation and corticospinal excitability has been investigated, the dose or energy of the stimulation has not attracted much attention and might play a pivotal role in modulating the activity of sensory regions. There is evidence showing that the number of peripheral stimulation pulses (i.e., dose) over time and the inter-pulse interval influence ongoing neural oscillations and have a neuromodulatory effect on the somatosensory cortex affecting indirectly corticospinal connectivity (Nitsche and Paulus, 2000;Schaworonkow et al., 2018;Zrenner et al., 2018). The nervous system maintains its excitability within an equilibrium range through adjustments derived from the history of neuronal activity, preventing excessive inhibition or facilitation (Pozo and Goda, 2010). Prolonged periods of stimulation-induced excitability/inhibition have been shown to activate homeostatic plasticity mechanisms that drive the system toward a more inhibited/excited state, reducing the effects of the stimulation (Gamboa et al., 2010;Andrews et al., 2013). Further, a progressive perceptual adaptation (or reduction of sensory responsiveness) has been evidenced after prolonged intervals of peripheral vibrotactile and electrocutaneous stimulation (Leung et al., 2005;Buma et al., 2007;Graczyk et al., 2018), indicating a neural compensation after a perturbation of the oscillatory neural system. Therefore, it could be conceivable that prolonged periods of NMES result in changes of neural excitability, resulting in a rebalance of the cortical response along time. Although most of the studies so far have used fMRI and NIRS to track slow brain correlates, or TMS-induced MEPs to asses corticospinal excitability, electroencephalography (EEG) is a powerful neuroimaging technique with high temporal resolution that has been used to study sensorimotor processes (including sensory evoked potentials and rhythms) (Birbaumer et al., 1990;Buzsaki, 2006;Shibasaki and Hallett, 2006;Leodori et al., 2019). The sensorimotor oscillations that mainly comprise the rolandic alpha [(7–13) Hz] and beta [(14–30) Hz] rhythms have been thoroughly used to study cortical involvement during sensorimotor tasks (Ramos-Murguialday and Birbaumer, 2015; López-Larraz et al., 2018), being quantified as the eventrelated (de)synchronization (ERD/ERS) (Pfurtscheller and Lopes da Silva, 1999). Furthermore, it has also been used as a feature for neuromodulation of sensorimotor neural network via proprioception and haptics (Ray et al., 2020;SebastiánRomagosa et al., 2020). To date, only a few studies have reported how oscillatory activity measured with EEG is modulated by NMES (Vidaurre et al., 2016, 2019;Tu-Chan et al., 2017;Corbet et al., 2018). Understanding this process is of great importance given the relevance of EEG and NMES for neurorehabilitation, especially for EEG-based neural interfaces. It is important to note that, although the neural response can be directly influenced by every stimulation parameter (e.g., frequency, pulse width, and pulse waveform), we focused on investigating the influence of intensity and dose of stimulation, since these are the parameters that are usually personalized for each patient. In this work, we acquired EEG activity from 12 healthy participants during NMES of the wrist extensor muscles at three different intensities (two below and one above the motor threshold) to investigate the neuromodulatory effect of peripheral stimulation on the ongoing cortical oscillatory activity, while the subjects rested in a comfortable sitting position. Firstly, we wanted to confirm that as the intensity increases, there is a proportional neural excitation, probably resulting from the recruitment of more afferent fibers, which provide a greater projection of afferent volleys to the sensory cortex. Secondly, we wanted to assess if the dose of stimulation has an effect on the magnitude of neural excitation over time, expecting that after the stimulation has been provided during a prolonged time, the cortical responsiveness is lower. Frontiers in Neuroscience | www.frontiersin.org 2January 2021 | Volume 14 | Article 593360 fnins-14-593360 January 15, 2021 Time: 14:6 # 3 Insausti-Delgado et al. NMES Parameters Influence Cortical Excitability MATERIALS AND METHODS Participants Twelve right-handed healthy participants (four females, age = 27.5 ±3.0 years, height = 176.3 ±8.6 cm, weight = 69.8 ±9.8 kg) were recruited to participate in the study. All of them signed an informed consent form. The experimental procedure was approved by the Ethics Committee of the Faculty of Medicine of the University of Tübingen (Germany). Participants were asked to stay comfortably seated on a chair with their right arm resting on a side table and the hand hanging with the palm facing downwards. Neuromuscular electrical stimulation (NMES) electrodes were placed on the righthand extensors (more detailed description in “Neuromuscular Electrical Stimulation” section), as described in Figure 1A. EEG and muscle activity of two sensors were recorded during the experiment. The electrical artifact recorded from the muscle sensors was used to align stimulation onset during the EEG signal processing. Experimental Design and Procedure The main purpose of the experiment was to investigate NMES neuromodulatory effects (instantaneous and cumulative) on brain oscillatory activity. With this aim, we compared the afferent cortical activity generated by three different NMES intensities. Participants were passively stimulated, meaning that they were resting, and no volitional motor command was generated during stimulation. Each participant underwent one session consisting of nine blocks, each comprising 18 trials. One of the three NMES intensities was randomly assigned to each block (determination of the current intensities explained in “Neuromuscular Electrical Stimulation” section), resulting in three blocks per intensity. A ready cue was presented 2.6– 3 s before the NMES interval, which had a random duration between 3.4 and 3.8 s. From the offset of the NMES to the next ready cue, a 3 s inter-trial period was introduced (see Figure 1B). Auditory cues announced the beginning of each interval. The time between blocks was used as breaks, lasting around 150 s (i.e., 2.5 min). The entire session including setup did not exceed 90 min. Data Acquisition The electroencephalographic activity was recorded with a commercial 32-channel actiCAP system (Brain Products GmbH, Germany) and a monopolar BrainAmp amplifier (Brain Products GmbH, Germany). The recording electrodes were placed at FP1, FP2, F7, F3, Fz, F4, F8, FC3, FC1, FCz, FC2, FC4, C5, C3, C1, Cz, C2, C4, C6, CP5, CP3, CP1, CPz, CP2, CP4, CP6, P7, P3, P4, P8, O1, and O2, following the international 10/20 system. Ground and reference electrodes were placed at AFz and Pz, respectively. Muscle activity of the right forearm of the participants was recorded by an MR-compatible BrainAmp amplifier (Brain Products GmbH, Germany) using two Ag/AgCl bipolar sensors (Myotronics-Noromed, Tukwila, Wa, United States). The sensors were placed laterally to the stimulation pads (see Figure 1A), using the right collarbone as ground. Both EEG and muscle activity were synchronously acquired at a sampling rate of 1,000 Hz. Neuromuscular Electrical Stimulation A programmable neuromuscular stimulator Bonestim (Tecnalia, Serbia) was used to deliver the stimulation. The cathode (3 ×3.5 cm, self-adhesive electrode) was placed over the muscles involved in wrist extension (extensor digitorum and extensor carpi ulnaris), at one third of the distance between the lateral epicondyle of the humerus and the Lister’s tubercle at the wrist. The anode (5 ×5 cm, self-adhesive electrode) was placed 5 cm distal to the cathode. To ensure the correct location of the electrodes, individually determined for each participant, stimulation above the motor threshold was applied until a complete wrist extension was induced. The frequency of the NMES was set to 35 Hz, and the pulse width to 300 µs (Lynch and Popovic, 2008). The individual intensities for each subject were obtained by a scan of currents, starting at 1 mA and increasing in steps of 1 mA. The participants were asked to report the initiation of the following sensations: (i) tingling of the forearm (i.e., sensory threshold—STh), (ii) twitching of the fingers (i.e., motor threshold—MTh), and (iii) complete extension of the wrist (i.e., functional threshold—FTh). According to these thresholds, the three NMES intensities (two below and one above the motor threshold) were calculated, following the equations defined in the study of Smith et al. (2003). Low intensity was defined as one third between the STh and MTh; medium intensity as two thirds between the STh and MTh; and high intensity as the FTh (see Eq. 1–3 and Figure 1C). These intensities were kept constant throughout the experiment, and the complete extension of the wrist induced by functional stimulation was visually verified. None of the participants reported pain or any harmful effect due to the stimulation. Low intensity =MTh −STh×0.33 +STh (1) Medium intensity =MTh −STh×0.66 +STh (2) High intensity =FTh (3) Data Preprocessing and Analysis Artifact Removal Procedures One important limitation for the quantification of EEG activity during continuous stimulation is the contamination of the signals due to the electrical currents delivered to the body. The EEG is easily polluted by these currents, and artifact removal methodologies are essential to properly estimate cortical activation. With this aim, different techniques for contamination removal in invasive and non-invasive brain activity recordings have been proposed, such as interpolation, blanking, or linear regression reference (LRR) (Walter et al., 2012;Iturrate et al., 2018;Young et al., 2018). Blanking of the data is the most restrictive method as contaminated data are rejected and signals that could be of interest are neglected for further analysis. Frontiers in Neuroscience | www.frontiersin.org 3January 2021 | Volume 14 | Article 593360 fnins-14-593360 January 15, 2021 Time: 14:6 # 4 Insausti-Delgado et al. NMES Parameters Influence Cortical Excitability Bonestim Intensity Sensory threshold Motor threshold Functional threshold 33% 66%0 Low intensity Medium intensity High intensity Start ReadyReady NMES 3 s2.6 ~ 3 s 0.8 s 3.4 ~ 3.8 s AB C FIGURE 1 | Experimental design and procedure. (A) Representation of the location of stimulation electrodes and sensors to measure muscle activity placed on right wrist extensors. (B) Timeline of the three phases included in each trial: preparation, NMES, and inter-trial period. (C) Determination of NMES intensities: low intensity as one third between the sensory threshold and motor threshold, medium intensity as two thirds between the sensory threshold and motor threshold, and high intensity as functional motor threshold. However, if the removal is implemented using hardware, the artifact has less influence on the recovery period of the amplifier, preventing it from being saturated and allowing the use of other methods to compensate for the missing data (Kent and Grill, 2012). Another approach is to linearly interpolate the corrupted data, connecting the last point before the artifact and the first point after the artifact. However, interpolation induces a bias in the estimation of power spectrum of the signals (Walter et al., 2012). LRR re-references the signals through weights that are assigned to each channel. The weights are calculated in a training block according to the noise of each channel generated by the electrical stimulation. This method effectively reduces artifacts, but it fixes the weights and assumes no changes in channel noise during the intervention. So far, the feasibility of this method has only been proven in invasive recordings (Young et al., 2018), in which impedances are less likely to change within sessions and are more similar among channels (Ball et al., 2009). Normally, impedances deteriorate and noiseinfluence increases throughout an EEG session, complicating the implementation of LRR in non-invasive recordings of brain activity. Therefore, we implemented an alternative two-step artifact removal method and demonstrated its feasibility. The raw EEG signals were pre-processed using custom-developed scripts in MATLAB (MathWorks, Natick, MA, United States). Channel removal based on power-line noise During an EEG session, particularly during setup and during periods between experimental blocks, special care is required to maintain EEG signal clean (i.e., raw data inspection and impedance check). However, our empirical experience shows that, sometimes, certain EEG electrodes present higher contamination due to the stimulation than others (see Figure 2A). These electrodes present broadband artifacts that impede further analyses even after applying the median filter preprocessing described below. We hypothesized that this effect might be due to degraded impedances, which occasionally deteriorate even when during the setup were set below 5 kOhm. Despite we did not store the impedances of each electrode to check this and discard the electrodes with high impedance, we ideated an automatized method to detect and discard them offline. This way, we could automatically eliminate contaminated channels without the human bias that would constitute a manual rejection. Having high impedance between the recording electrode and the skin resembles an open circuit, where the electrode behaves like an antenna and captures outside electric frequencies (like the power-line noise). Our method exploits this effect and identifies the EEG electrodes with unusually high power-line noise (50 Hz in Europe). Since all the electrodes should be approximately equally exposed to electromagnetic signals at 50 Hz, we assume that very high power at this frequency is an indirect indicator of high skin-electrode impedance. The procedure was applied block-wise, meaning that it was used to detect and remove contaminated channels within each individual EEG block. The EEG activity was high-pass filtered at 0.1 Hz with a 4th-order Butterworth. The power spectral distribution at 48–52 Hz was estimated using Welch’s method, averaging the periodogram of 1 s Hamming windows with 50% overlapping. The power mean and standard deviation (SD) of all the EEG channels were calculated from all trials in each block. Channels whose power was higher than 4 SD above the mean were discarded from that specific block. The remaining channels were used to re-compute the mean and SD. The procedure was iteratively repeated until no channels exceeded the rejection threshold (see Figure 3 dashed box, Supplementary Figures 1,2, and Supplementary Table 1). Median filtering for removal of electrical stimulation contamination It is well known that applying electro-magnetic currents to stimulate the neural system can introduce undesired noise to the recordings. The NMES configuration used in this study introduces large peaks of short latency (∼5 ms) to the recorded EEG signal. Therefore, median filtering was used to minimize the NMES-induced artifacts (Insausti-Delgado et al., 2017). This filter is suited to eliminate high-amplitude peaks from a time series (Gallagher and Wise, 1981) and can remove the shortlatency high-amplitude artifacts caused by the NMES. A sliding window of 10 ms was applied to the EEG signal in steps of one Frontiers in Neuroscience | www.frontiersin.org 4January 2021 | Volume 14 | Article 593360 fnins-14-593360 January 15, 2021 Time: 14:6 # 5 Insausti-Delgado et al. NMES Parameters Influence Cortical Excitability 0123-1-2-3 44Time (s) C3 C1 Cz CP3 CP1 120 µV Contamination of channels with bad impedances Time (s) 2.10 2.12 2.14 2.16 2.18 2.20 C1 NMES artifact removal by median filter Raw signal Median filtered 15 µV AB FIGURE 2 | Characterization of contamination induced by NMES. (A) EEG of a representative trial showing non-contaminated channels (C3, C1, CP3, and CP1) and one contaminated channel (Cz) during high-intensity stimulation. Channels with bad impedances are more prone to be contaminated by electrical stimulation. (B) Zoom in 100-ms segment of a representative EEG trial in a non-contaminated channel that presents NMES artifacts (red line). The effect of stimulation artifacts is minimized by median filter (green line). Time-frequency analysis HP filter 0.1 Hz Power of EEG channels [48-52] Hz Compute mean and SD power over all channels Power EEG Ch > Mean power + 4SD ? EEG channels with good impedances Discard EEG channel YES NO Quantification of brain oscillatory activity Estimation of power in alpha (7-13) Hz and beta ( 14 -30) Hz Non-stimulation (-3, -1) s Low NMES (0.5, 2.5) s Med. NMES (0.5, 2.5) s High NMES (0.5, 2.5) s Raw data Channel removal Median filter BP filter 0.1-45 Hz Trial extraction (-4, 4) s Common average re-referencing Subsampling to 100 Hz FIGURE 3 | Flowchart with the steps for the quantification of brain oscillatory activity. The whole set of data is firstly preprocessed by a two-step procedure based on channel removal and median filtering. In this level, channels with bad impedances and artifacts due to electrical stimulation are removed. Then, the remaining clean data are filtered and divided into trials. Finally, power is estimated in alpha (7–13) Hz and beta (14–30) Hz bands for non-stimulation [(−3, −1) s] and NMES [(0.5, 2.5) s] intervals, being the baseline (−2.5, −1.5) s. Frontiers in Neuroscience | www.frontiersin.org 5January 2021 | Volume 14 | Article 593360 fnins-14-593360 January 15, 2021 Time: 14:6 # 6 Insausti-Delgado et al. NMES Parameters Influence Cortical Excitability sample, providing as output the median value of each window. We selected a 10 ms window as it fully covers the electrical artifact. This filter produces a frequency-dependent attenuation that follows an exponential function from 0 to 100 Hz (i.e., the frequency with a period that completely fits within the 10 ms window), leading to low attenuation at low frequencies and a complete attenuation at 100 Hz (Supplementary Figure 3). With this window size, the attenuation of the signal at 10 Hz, 20 Hz, and 30 Hz is 1.28%, 4.89%, and 10.90% respectively. The relatively low attenuation at low frequencies makes this method suitable for analyzing alpha and beta sensorimotor oscillations. Figure 2B displays a zoomed segment of 100 ms of activity, showing the effect of the median filter on the stimulation artifacts. Quantification of Brain Oscillatory Activity A common average reference (CAR) was applied to the EEG signals. The re-referenced signals were band-pass filtered at 0.1– 45 Hz using a 1st-order Butterworth filter. Each block was trimmed down to (18×) 8-s trials, from -4 to + 4 s, being 0 the beginning of the stimulation. The trials were down sampled at 100 Hz, and those belonging to the same level of NMES intensity were pooled together. The quantification of cortical activity was performed by evaluating the spectrum differences of the sensorimotor rhythms, by means of the alpha and beta event-related (de)synchronization (ERD/ERS), i.e., decrease or increase in power generated by an event compared to a baseline (Pfurtscheller and Lopes da Silva, 1999). Large ERD values (i.e., more negative power values) represent stronger cortical activation compared to baseline time interval, as it represents disinhibition/excitation of neural population activity (Ritter et al., 2009). For the quantification of brain activity, we used the FieldTrip toolbox1for MATLAB. Time-frequency maps were calculated using Morlet wavelets in the frequency range from 1 to 45 Hz, with a resolution of 0.5 Hz. The power change was computed as the percentage of increase or decrease in power (i.e., ERS or ERD) with respect to the baseline [(−2.5, −1.5) s], where Pjrepresents the signal power at the jth sample, as described in Eq. 4. ERD/ERSj(%)=Pj−Baseline Baseline ×100 (4) From the time-frequency maps, we calculated the sensorimotor averaged changes in power in alpha (7–13) Hz and beta (14– 30) Hz bands for non-stimulation [(−3, −1) s] and NMES [(0.5, 2.5) s] intervals, using as baseline the (−2.5, −1.5) s interval (see Figure 3). The NMES period was defined as starting at 0.5 s to avoid potential bias and influence of the stimulation onset like event-related brain potentials (e.g., somatosensory evoked potential) at t= 0. For the topographical inspection and the descriptive analysis of the entire brain activity, all EEG channels were analyzed individually. For the quantitative analysis, we calculated the mean change in power of channels C1, C3, CP1, and CP3 (i.e., area over the sensorimotor cortex 1http://fieldtriptoolbox.org/ representing the right forearm, contralateral hemisphere to the stimulated limb), since we considered that averaged values over these electrodes could better quantify the overall changes in the sensorimotor areas. Statistical Analysis The statistical tests were performed in IBM SPSS 25.0 Statistics software (SPSS Inc., Chicago, IL, United States) and MATLAB. We used the Shapiro–Wilk test to determine the normality of the data. Accordingly, a multivariate analysis of variance (MANOVA) for repeated measures was performed to find differences in the dependent variables, alpha and beta ERD/ERS, with NMES intensity (four levels: no stimulation, low-, medium-, and high-intensity stimulation) as within-subject factor. In order to determine the origin of the significant effect, post hoc tests with Bonferroni correction were performed. In order to analyze whether NMES can induce a doseeffect, we studied the ERD/ERS changes over time. For that, we computed the alpha and beta ERD/ERS for each single trial (i.e., in Eq. 4, Pjwas the alpha/beta power of each trial during the NMES period, and the baseline was calculated from the grand average of all the trials of each intensity). A linear regression was estimated for the ERD/ERS values over trials for the two frequency bands (i.e., alpha and beta) and the three NMES intensities (i.e., low, medium, and high). Correlation between ERD/ERS and sequence of trials were calculated using Pearson’s correlation coefficient to study stimulation effects over time. RESULTS Effect of Artifact Removal The pre-processing of the data eliminated satisfactorily the electrical noise contamination coming from the peripheral electrical stimulation. It reduced the effect of the artifacts to an extent that allowed us to perform EEG spectral analysis of the brain oscillatory activity. Channels with good impedances are also influenced by the electrical stimulation artifact that is introduced into the signal as large peaks. Figure 2B illustrates in a 100 ms segment of a representative trial how the median filter deals with these undesired artifacts. To prove the efficacy of the method, we focused on the worst-case scenario, as the contamination is larger for higher stimulation intensities. This effect can be observed in Figure 4, which depicts the EEG time-frequency activity at the different NMES intensities including artifacts and after the median and spatial filters are applied. The NMES generates an increase of power, or ERS, around 35 Hz (i.e., the stimulation frequency), which increases as the NMES intensity is incremented (Figure 4, left column). This power increase in highbeta/low-gamma band due to stimulation artifact was eliminated for all intensities after median filtering. Applying the median filter did not change the power in alpha and beta frequencies before the stimulation onset (t= 0 s), but eliminated the ERS during the stimulation period, minimizing the artifacts and revealing the alpha and beta modulation. The common average re-reference (CAR) after median filtering enhanced the power decrease of Frontiers in Neuroscience | www.frontiersin.org 6January 2021 | Volume 14 | Article 593360 fnins-14-593360 January 15, 2021 Time: 14:6 # 7 Insausti-Delgado et al. NMES Parameters Influence Cortical Excitability Without median filter With median filter 5 15 25 35 45 0123-1-2-3 Time (s) Frequency (Hz) 5 15 25 35 45 0123-1-2-3 Time (s) Frequency (Hz) 50400102030-50 -40 -30 -20 -10 ERD/ERS (%) 5 15 25 35 45 0123-1-2-3 Time (s) Frequency (Hz) 5 15 25 35 45 0123-1-2-3 Time (s) Frequency (Hz) 5 15 25 35 45 0123-1-2-3 Time (s) Frequency (Hz) 5 15 25 35 45 0123-1-2-3 Time (s) Frequency (Hz) Medium intensity Low intensity High intensity 5 15 25 35 45 0123-1-2-3 Time (s) Frequency (Hz) 5 15 25 35 45 0123-1-2-3 Time (s) Frequency (Hz) 5 15 25 35 45 0123-1-2-3 Time (s) Frequency (Hz) With median filter and CAR FIGURE 4 | Comparison of cortical activation after median and spatial filtering. Time-frequency maps averaged over all participants, representing ERD/ERS, of the average of channels (C3, C1, CP3, and CP1) located over the contralateral sensorimotor cortex to the stimulated limb. Averaged time-frequency maps without median filter (left column), with median filter (center column), and with median and CAR filter (right column) after removal of contaminated channels. Rows show the different NMES intensities: high (upper row), medium (middle row), and low intensity (lower row). The percentage of ERD/ERS is computed according to the baseline (−2.5, −1.5) s. Time 0 s is aligned with the onset of the stimulation. the bands of interest for every NMES intensity, as it does with non-contaminated EEG (Wolpaw et al., 2002). Regardless of the intensity delivered, the initiation of the stimulation generated timeand phase-locked activity (i.e., event-related potentials— ERP), presumably due to the sensory processing of the NMES (Reynolds et al., 2015). This can be seen as power increase at low frequencies (1–5 Hz). Influence of Stimulation Intensity on Cortical Activation To analyze the influence of stimulation intensity on cortical activation, we compared the changes in brain oscillatory activity in four conditions: non-stimulation, low-, medium-, and high-intensity stimulation (see topoplots of alpha and beta rhythms in Figure 5). Topographic maps of non-stimulation condition were calculated using the interval (−3, −1) s prior to the stimulation, while the other conditions were extracted from the interval (0.5, 2.5) s after stimulation onset. An increment of the stimulation intensity resulted in an increasing ERD (i.e., larger decrease in power) in both frequency bands over the sensorimotor cortex as expected (Backes et al., 2000;Smith et al., 2003;Schürholz et al., 2012), while occipital areas showed idling activity. At high-intensity stimulation, the sensorimotor cortex of both hemispheres presented a decrease of power, being more pronounced in the contralateral hemisphere, as demonstrated in previous work studying brain oscillatory signatures of motor tasks (Ramos-Murguialday and Birbaumer, 2015). Frontiers in Neuroscience | www.frontiersin.org 7January 2021 | Volume 14 | Article 593360 fnins-14-593360 January 15, 2021 Time: 14:6 # 8 Insausti-Delgado et al. NMES Parameters Influence Cortical Excitability Medium intensity ytisnetnihgiHytisnetniwoLNon-stimulation Alpha band (7-13 Hz) 40020-40 -20 ERD/ERS (%) * ** 0 -10 -20 ERD/ERS (%) -30 NonStim. Low Med. High 0 -10 -20 ERD/ERS (%) -30 NonStim. Low Med. High * ** Beta band (14-30 Hz) FIGURE 5 | Comparison of cortical activation for different stimulation intensities for alpha and beta band EEG activity. Topographic maps, averaged over all participants, showing ERD/ERS of non-stimulation periods (−3, −1) s and NMES periods (0.5, 2.5) s belonging to each intensity (i.e., low, medium, and high) for alpha (upper row) and beta (lower row) frequency bands. Bar graphs show the mean percentage of ERD/ERS averaged from channels (C3, C1, CP3, and CP1) for each intensity and frequency band. The statistically significant differences between pairs are expressed with horizontal lines and stars. The percentage of ERD/ERS is calculated with respect to the baseline (−2.5, −1.5) s. The signals are processed using the two-step procedure (i.e., removal of contaminated channels and median filter) and a CAR. Our MANOVA analysis reflected a significant effect of intensity on alpha and beta ERD [F(6, 86) = 4.356, p= 0.001]. Rightmost panels in Figure 5 display the results of the post hoc comparisons. For both alpha and beta, there was a significantly higher ERD (i.e., more negative power values) induced by high-intensity NMES compared to rest (p= 0.004 for alpha, p<0.001 for beta), to low-intensity NMES (p= 0.013 for alpha, p= 0.004 for beta), and to medium intensity (p= 0.045 for alpha, p<0.001 for beta). Stimulation Dose-Effect To study the influence of stimulation dose on cortical activity, we computed the ERD/ERS of each single trial and performed a regression in time within session appending same stimulation intensity blocks in order of appearance (blocks of different intensities were presented randomly). Figure 6 shows the average ERD/ERS for all participants in both frequency bands during the 54 trials (3 blocks ×18 trials, see vertical yellow lines) for each intensity and a linear regression to fit them. We also performed regressions over time within each block with consecutive trials (see Supplementary Figure 4). The first thing we observed is a clear modulation based on stimulation, reducing its variability with stimulation amplitude. Low and medium intensities caused a significant reduction of alpha ERD over time (p= 7.8e-4 for low; p= 0.0053 for medium). In contrast, high intensity caused a significant enhancement of beta ERD over time (p= 5.12e-5). Furthermore, as can be seen in Figure 6, we observed that in every block (separated by yellow vertical lines), there is a reduction of the ERD (increase of power plotted as linear regression in Supplementary Figure 4) progressively induced per trial from the first to the last trial. The first trial of the new block presented a larger ERD (decrease of power) in comparison to the ERD of the last trial of the previous block (irrespective of the stimulation intensity and order of the block within the session). DISCUSSION This study demonstrated the significant effects of artifacts, intensity, and dose of neuromuscular electrical stimulation (NMES) of the upper limb on the ongoing brain oscillatory activity recorded using EEG. First of all, we dealt with the issue of artifact removal to allow accurately estimating the cortical oscillatory activity. Recordings of brain activity are easily polluted, especially when electrical stimulation interacting with the nervous system is concurrently used. This contamination can negatively affect the signal-to-noise ratio, covering the brain Frontiers in Neuroscience | www.frontiersin.org 8January 2021 | Volume 14 | Article 593360 fnins-14-593360 January 15, 2021 Time: 14:6 # 9 Insausti-Delgado et al. NMES Parameters Influence Cortical Excitability Low intensity 0 20 -20 40 -40 ERD/ERS (%) 0504030102 Number of trial ytisnetnihgiHytisnetnimuideM Alpha band (7-13 Hz) 0 20 -20 40 -40 ERD/ERS (%) 0504030102 Number of trial 0 20 -20 40 -40 ERD/ERS (%) 0504030102 Number of trial 0 20 -20 40 -40 ERD/ERS (%) 0504030102 Number of trial 0 20 -20 40 -40 ERD/ERS (%) 0504030102 Number of trial 0 20 -20 40 -40 ERD/ERS (%) 0504030102 Number of trial Beta band (14-30 Hz) R = -0.5222, m = -0.2374, p = 5.12e-05R = -0.1305, m = -0.0687, p = 0.3469R = 0.2189, m = 0.1025, p = 0.1116 R = 0.4433, m = 0.4564, p = 7.8e-04 R = 0.3740, m = 0.3294, p = 0.0053 R = 0.0931, m = 0.0610, p = 0.5030 FIGURE 6 | Comparison of cortical activation over trials for alpha and beta bands. The cortical activity during NMES period [(0.5, 2.5) s] quantified as ERD/ERS over the 54 trials for each stimulation intensity, divided into blocks by vertical yellow lines, and averaged for all the participants. The percentage of ERD/ERS is calculated according to the baseline (−2.5, −1.5) s. Different intensities are compared in columns: low (left), medium (middle), and high (right). Alpha (upper row) and beta (lower row) frequency bands are described. Significant correlations between ERD/ERS and sequence of trials over session are represented with black solid linear regressions. activity. Our findings evidenced that the median filter enhanced the detection of sensorimotor oscillatory activity after removing stimulation artifacts. After the EEG data was cleaned, especially of NMESinduced artifacts, we analyzed the modulation of alpha and beta oscillations produced by the stimulation. Power suppression, or desynchronization, of these frequencies has been associated with cortical excitation, whereas synchronization reflects a state of inhibition (Klimesch et al., 2007). During highintensity NMES, the induced desynchronization in alpha and beta was significantly larger than during stimulation at low or medium intensities or no stimulation. While below-motor-threshold stimulation intensities only activate cutaneous mechanoreceptors (e.g., Pacinian corpuscles and Merkel disks) and sensory axons, stimulation above the motor threshold also recruits proprioceptive receptors (e.g., muscle spindles, Golgi tendon organs, and joint afferents) (Maffiuletti et al., 2008;Bergquist et al., 2011;Golaszewski et al., 2012). It has been proposed that muscle spindles transmit inputs to the spinal cord and can directly influence the motor cortex (M1) through the area 3a, while the projections from area 3b activated by cutaneous feedback to M1 are less likely to happen (Carson and Buick, 2020). We can therefore assume that high-intensity NMES leads to higher neural excitation, probably by recruiting a larger number of receptors derived from muscle contractions, in addition to the cutaneous and sensory fiber afference that is also engaged in sensory-threshold stimulation, and that the recruitment of muscle spindles results in activation of M1 via area 3a of the somatosensory cortex (S1) (Schabrun et al., 2012;Carson and Buick, 2020). These results of intensitydependent brain activation are in line with corticomuscular responses (Sasaki et al., 2017), metabolic responses recorded by functional magnetic resonance imaging (fMRI) (Backes et al., 2000;Smith et al., 2003) and near-infrared spectroscopy (NIRS) (Schürholz et al., 2012), which demonstrated a direct quantitative association with stimulation intensity. Noteworthy, our results showed that the cortical activity measured with EEG, quantified as event-related (de)synchronization (ERD/ERS), during low and medium intensities was not significantly different to no stimulation. This suggests that below-motorthreshold NMES might not recruit enough afferent fibers to induce significant cortical modulation and that more afference (probably through muscle contraction, proprioception due to the movement, and a larger number of sensory fibers recruited) is required to transmit more information that reaches the brain and is measurable in the EEG at the analyzed frequencies. It is well known that during voluntary movement, a stronger cortical activation is seen in alpha than in beta (López-Larraz et al., 2014;Ramos-Murguialday and Birbaumer, 2015). Such modulation in alpha and beta cortical activities has been related to the control of top-down and bottomup neural processes, suggesting its role in the integration Frontiers in Neuroscience | www.frontiersin.org 9January 2021 | Volume 14 | Article 593360