Gamma-band enhancement of functional brain connectivity following transcutaneous electrical nerve stimulation

Objective. Transcutaneous electrical nerve stimulation (TENS) has been suggested as a possible non-invasive pain treatment. However, the underlying mechanism of the analgesic effect of TENS and how brain network functional connectivity (FC) is affected following the use of TENS is not yet fully understood. The purpose of this study was to investigate the effect of high-frequency TENS on the alteration of functional brain network connectivity and the corresponding topographical changes, besides perceived sensations. Approach. Forty healthy subjects participated in this study. Electroencephalography (EEG) data and sensory profiles were recorded before and up to an hour following high-frequency TENS (100 Hz) in sham and intervention groups. Brain source activity from EEG data was estimated using the LORETA algorithm. In order to generate the functional brain connectivity network, the Phase Lag Index was calculated for all pair-wise connections of eight selected brain areas over six different frequency bands (i.e. δ, θ, α, β, γ, and 0.5–90 Hz). Main results. The results suggested that the FC between the primary somatosensory cortex (SI) and the anterior cingulate cortex, in addition to FC between SI and the medial prefrontal cortex, were significantly increased in the gamma-band, following the TENS intervention. Additionally, using graph theory, several significant changes were observed in global and local characteristics of functional brain connectivity in gamma-band. Significance. Our observations in this paper open a neuropsychological window of understanding the underlying mechanism of TENS and the corresponding changes in functional brain connectivity, simultaneously with alteration in sensory perception.


Introduction
Pain and sensation are complex phenomena and the mechanism of how the brain process the sensory and painful input is not well understood. It has been shown that five main cortical regions play a significant role in pain and sensory regulation, namely: the primary somatosensory cortex (SI), insular cortex, the anterior cingulate cortex (ACC), the secondary somatosensory cortex (SII), and the medial prefrontal cortex (mPFC) [1][2][3].
Moreover, the oscillatory activity of brain areas has been highlighted in the somatosensory cortex in pain and sensory processing [4]. For example, magnetoencephalography (MEG) studies suggested dynamic brain activity in alpha, beta, and gamma bands as a signature for perception [5]. Oscillation in high-frequency electroencephalography (EEG) (i.e. beta and gamma bands) are speculated to represent bottom-up, feed-forward processing that is sensorydriven. However, lower frequency bands synchronization are believed to represent regulation of brain activity as well as top-down feedback processing [6].
Transcutaneous electrical nerve stimulation (TENS) is a neurorehabilitation and neuromodulation technique for the induction of pain relief as it is believed to affect the nervous system both at the peripheral and central levels [7][8][9]. The underlying mechanism of TENS at the central level has been reported by alteration in the sensory system [10].
Sensory and pain-related evoked potentials have shown alterations of the sensory responses in the central nervous system, caused by TENS, through inhibition of the SI activity [10]. The lasting effect of TENS on the somatosensory cortex has been investigated recently [11], and results demonstrated a significant decrease in the N100 and P200 waves' amplitudes, at least up to one hour following TENS intervention. In addition, a statistically meaningful suppression in the perceived sensation intensity as the effect of TENS has been reported. Moreover, studies have been shown that the analgesic effect of TENS depends on different waveform properties of the stimulation pulses (including frequency, intensity, and pulse width) [12][13][14]. For example, the application of TENS with the intensity above the motor threshold resulted in an increase in corticomotor excitability. The corticomotor excitability was suppressed when an intensity below the motor threshold was used [15,16].
The graph theory combined with advanced biosignal processing has been proven as a powerful tool to describe the functional synergistic behavior of various brain regions for different conditions [17][18][19]. In this regard, the activity of brain regions is considered as the nodes of the network, and the strength of connectivity between two nodes defines the vertices. To quantify the strength of communication, amplitude-based and phase-based indices have been used based on data collected by several neuroimaging techniques (e.g. fMRI [20,21], EEG [19,22], MEG [23,24]). The quality and intensity of the information exchange within a network can be mathematically explained using graph theory. Previous studies have suggested 'small-world topological features (such as clustering coefficient, degree, global efficiency)' of the network as strong indicators and predictors of brain alteration in several neurological disease patients (please see [25] and references therein).
The effect of TENS on pain and sensation may be caused by the alteration in the brain network, and the corresponding functional connectivity (FC) features between various brain regions [26]. In this paper, motivated by the lack of understanding of the underlying neurophysiological mechanism of TENS in pain therapy, we hypothesize that short-term and longterm effects of TENS on sensation can be decoded by analyzing the topographical alterations in the FC of the central nervous system. In this work, we study a group of 40 healthy subjects following non-noxious sensory input. Figure 1 illustrates the executed protocol. Four somatosensory evoked potentials (SEPs) were included in each experiment session. This was planned to evaluate the cortical changes due to TENS intervention in the following phases: Pre-SEP phase (T0)-considered as baseline. Three post-intervention phases were conducted immediately after intervention (T1), 30 min (T2), and 60 min (T3). At each of the four recording phases (T0, T1, T2, and T3) and for both groups (TENS and sham), we applied 40 double-pulse stimuli twice (a total of 80 trials for each SEP phase).

Participants
Forty healthy participants (20 female), 26.9 ± 4.3 [mean ± std] years old, were randomly recruited and assigned equally to an intervention group (n = 20, 10 female) or a sham group (n = 20, 10 female). Each participant in the study signed an informed consent form in accordance with the Committee on Health Research Ethics (N-20 180 049) in the North Denmark Region.

Data collection
During all experimental conditions, participants were seated in a comfortable chair in a temperaturecontrolled room (i.e. temperature between 24 • C and 26 • C). They were asked to focus on the stimuli and gaze at a black cross in the center of a screen. A 64-channel system was used using the international 10-20 system to record continuous EEG data. The sampled data were amplified and recorded using the BrainAmp MRI-compatible amplifier (Brain Products GmbH, Munich, Germany). All electrodes were grounded to an electrode that was placed between the Fz and Fpz electrodes and referenced to the FCz electrode. EEG was digitalized at a 5 kHz sampling rate (electrode impedances maintained below 20 kΩ).
Two succeeding constant-current square-wave electrical pulses with an inter-pulse interval of 10 ms (500 µs pulse width) were delivered to the left median nerve for all SEP phases [27]. Electrical stimulations were applied with an inter-stimulation interval randomly changed between 6 and 8 s to avoid habituation [28]. These electrical pulses were generated using a constant current stimulator (DS5, Digitimer, UK) and applied to two surface electrodes (Axelgaard PALS Electrodes, skin contact size 4 × 4.6 cm, oval) placed on the left-median nerve close to the wrist. The placement of the electrodes is given in figure 1.
As defined by the staircase procedure, the stimulus intensity for all SEP phases was individually adjusted to an amplitude that was twice the sensation threshold (without muscle contraction) [29]. A 0.5 mA initial intensity was considered for the electrical stimuli. This amplitude was iteratively increased by 0.5 mA until the stimulus was reportedly perceived by the subject. As soon as the subject perceived the stimulus, the stimulus intensity was reduced in steps of 0.3 mA until no sensation was perceived. Then, the increasing steps of current intensity adjusted to Figure 1. Experimental procedures. The brain network alteration following TENS was assessed by analyzing EEG data recorded at four time points (T0 as baseline recording, T1 as immediately after the intervention, T2 as 30 min, and T3 as 60 min after intervention). TENS and sham interventions was applied following T0 (baseline recording) for each group. Surface electrodes (Axelgaard PALS Electrodes, skin contact size 4 × 4.6 cm, oval) placed on the left-median nerve close to the wrist. 0.1 mA and lasted until the stimulus was reported again by the subject. The sensation threshold was defined as the average intensity of the measured intensities through three times staircase procedure. This procedure was used for participants subjected to both TENS and Sham intervention.
Subjects were asked to rate the sensation intensities perceived by sensory electrical stimulation in the SEP phases using a 0-10 numerical rating scale from no touch to maximum non-painful sensation. Additionally, following a block of SEP recording, individual perception profiles were recorded by asking subjects to highlight the area of the perceived sensation with the relevant intensity on custom-made software on the screen.

TENS intervention
We applied a high-frequency TENS pattern (i.e. defined as a 100 Hz pulse train with a high-intensity amplitude that however does not cause motor movement) similar to what has been reported by Lai et al and Cruccu et al in their respective studies pertaining to acute and chronic pain, and stroke rehabilitation [7,9]. The electrical stimulation for both TENS and Sham intervention was applied for 20 min. For TENS, it includes 40 repetitions of the following phases: a 20 s on-stimulation phase (during which 100 Hz trains of 1 ms pulse width were provided) and a 10 s off stimulation interval. For the intervention group, without causing motor response or discomfort, the stimulation intensity was set at 80% of the discomfort level. The stimulation intensity was adjusted to be the same as the sensation threshold for subjects in the sham-controlled group. For the Sham group, the stimulation lasted for 1 min (i.e. two trials of stimulation), and for the remaining 19 min, high-frequency electrical stimulation was not delivered to the subjects. The participants were informed: 'For the next 20 min, the electrical stimulation would deliver to your median nerve. The perceived sensation will be subjective and range from no or weak sensations to intense sensations' .

Data pre-processing
Off-line EEG data were analyzed by BrainVision Analyzer 2.2 software (Brain Products® GMBH) and Matlab. Raw EEG data were downsampled to 256 Hz and re-referenced to the averaged reference. Next, data were band-pass filtered from 0.5 to 90 Hz (4th order Butterworth Zero Phase IIR Filter), and band-rejected filtered using a 2nd order Butterworth notch filter at 50 ± 1 Hz. EEG filtered data were then segmented into epochs of 2000 ms for each stimulation (including 500 ms before and 1500 ms after the stimulus onset). Blink correction of EEG signals using ICA has been reported to affect the amplitude and phase of the EEG signal [30][31][32]. Thus, in this work ICA-based correction was not conducted; instead, epochs contaminated by blink artifact, flagged using ICA analysis, and thus, on average, 15 epochs (out of total 80 epochs per subject) were deleted as part of signal inspection before advanced processing. Artificial power effects and phase synchrony activated by micro-saccadic eye movements were corrected using Laplacian current source density (CSD). CSD also decreased the volume conduction effect on phase synchronization [33]. Baseline correction was performed for individual epochs for the period of 500 ms before stimulus onset. Next, data were inspected to remove any remaining non-physiological artifacts and epochs with amplitude values exceeding ±100 µv or contaminated by excessive noise (3 ± 2 out of on average 65 epochs, following blink rejection) were excluded.
Individual SEP was measured by averaging the preprocessed epochs for each time phase. To confirm the validity of the applied methodology for the SEPs, the grand average of the individual SEPs for different time phases (T0, T1, T2, and T2) and conditions (TENS and Sham) are extracted at the Cz channel and illustrated in figure 2. Two SEP components, N2 and P2 waves, are represented in figure 2 as the most negative and positive peaks between 100 and 400 ms after stimulus onset, respectively.

Source localization
The source-level electrophysiological activity was estimated using the build-in low-resolution electromagnetic tomography (LORETA) algorithm of Brainvision analyzer 2.2. The LORETA method has been widely used to investigate brain activity [34]. The LORETA algorithm estimates the electrical neuronal activity in source space (i.e. current density values [µV 2 mm 4 ]) by computing an inverse solution. As such, 2394 voxels of spectral densities with 7 mm spatial resolution were computed using the MNI-305 brain template from a 64-channel scalp EEG [35].
Following anatomical regions of interest (ROIs) were specified using automated anatomical labeling (AAL) brain template [36] and specified in line with the sensory/pain processing network hypothesis described in the introduction. We selected eight ROIs, including the primary sensory cortex (SI_L and SI_R, on the left and right hemispheres, Precentral gyrus in AAL), the secondary sensory cortex (SII_L and SII_R, on the left and right hemispheres, Rolandic operculum in AAL), the mPFC (Superior frontal gyrus medial in AAL), the ACC, and the anterior insula (INS_L and INS_R, on the left and right hemispheres, Insula in AAL) and performed LORETA algorithm. The average current density value across all the voxels within an ROI resulted in a time series defined as the source activity.

Functional connectivity
All the data from selected sources were first filtered into six predefined frequency bands (δ: 0.5-4 Hz θ: 4-8 Hz, α: 8-13 Hz, β: 14-40 Hz, γ: 51-90 Hz and 0.5-90 Hz [37]) using a 2nd order Butterworth filter. The phase from a time-series signal was obtained by the Hilbert transform. Next, the FC was calculated for pair-wise connections between ROIs activities. The phase lag index (PLI) was calculated for each pair of ROI across segments for each subject, frequency bands, and time phases.
The PLI has been widely used to calculate FC in both the channel and source domain. As a reflection of the phase synchronization, PLI measures the degree of inter-trial phase variability [38]. Following disregarding and attenuating the zero and π phase differences, PLI suppresses the spurious synchronization from the volume conductance and the common sources of signal interference [39]. Connectivity between signals was therefore indicated by a divergence of symmetric distribution. The PLI was calculated as follows [40]: where ∆ϕ is the phase difference of two-time series (source activity) at a given time point (t k ). The extracted PLIs were then averaged across trials, and a synchronization matrix was created for each participant, time phase, and frequency band of all defined connections.

Graph analysis
The induced changes in the brain network following TENS intervention were also examined by graph theory analysis. Here, the eight selected ROIs described the nodes. The extracted PLI indexes defined the edges between the nodes, which resulted in a weighted and undirected network. For each constructed PLI network at each frequency band, the network alternations following TENS and Sham intervention were examined by node strength, clustering coefficient, and efficiency using the Brain Connectivity Toolbox [41]. Global efficiency represents the level of functional integration in the brain network, and this measure has been used in several studies to assess the changes in brain network FC [19,42]. Moreover, the extent of the network's functional segregation is explained by the global clustering coefficient [18,43], while the global strength depicts the overall network connection strength. Node strength is determined as the sum of the edge weights connecting the selected node and represents the strength of interconnectivity with other nodes [44]. The clustering coefficient is a metric indicating functional discrimination and is calculated by the fraction of the node's neighbors that are also connected to each other. The efficiency is described as the average inverse shortest path length between nodes, which indicates the performance of information exchange in a network.

Statistical analysis
The statistics were divided into two groups and conducted in RStudio (R version 4.0.3 and RStudio version 1.3.1093). First, statistical tests were applied to identify FC changes (i.e. for each edge) across time (i.e. T0, T1, T2, and T3) and group. Since Shapiro-Wilk test results (examined the normality of distribution) revealed that data in the first group did not follow a normal distribution, non-parametric tests were performed. In this regard, separately for each group (TENS and Sham), Friedman tests were employed to compare the PLI changes of each vertex across four time phases and for four different frequency bands. In the case where a statistical significance was found, Mann-Whitney tests were used as a post hoc analysis to investigate the lasting effect of the intervention between TENS and Sham. For each frequency and condition, Bonferroni correction was performed to address the multiple comparisons problem. Second, the alteration in the constructed FC networks (PLI networks) was explored by analyzing the network parameters such as clustering coefficient, node strength, and network efficiency. In addition to the above, for the second type of analysis, the normality of data distribution of the network characteristics passed Shapiro-Wilk test. As a result, an independent sample t-test was performed to analyze the significance of the difference (changes in network parameters in T1 compared to the T0) in the TENS and Sham groups. Moreover, bootstrap independent-sample t-tests with 1000 repetitions were applied to address robustness and multiple comparison problems. The level of significance was p < 0.05 in all tests.
Since Shapiro-Wilk tests revealed that perceived sensation data were not normally distributed, nonparametric statistical evaluations were used for data analysis. Friedman tests were employed to analyze the effect of TENS on sensation in four time phases (within-subject). Mann-Whitney tests were applied to compare perceived sensation data between TENS and Sham groups as between-subject factors. Multiple comparisons were addressed using Bonferroni correction.

Perceived sensation
The average sensation threshold at T0 in the TENS and Sham group was 2.39 ± 0.59 mA and 2.32 ± 0.59 mA, respectively. Statistical analysis was applied on the normalized perceived sensation (change in sensation compared to the T0 phase).
Results from the Friedman test on normalized perceived sensation revealed a statically significant difference between TESN (χ 2 (3) = 54.20, p < 0.001) and Sham (χ 2 (3) = 36.95, p < 0.001). Changes in sensation in T1, T2, and T3 following TENS application compared to T0 were statistically assessed with the Mann-Whitney test for the TENS and Sham group. Significant suppression in perceived sensation in the TENS group compared to the Sham group was found and remained up to an hour following TENS (T1-T0 p = 0.005, T2-T0 p = 0.002, T3-T0 p = 0.002).
The individual map of sensation profile following intervention was normalized to the recorded sensation profile at T0, and the normalized group average of hand sensation profile map was demonstrated to supplement the perceived intensity results (figure 3). Results show that, while the sensation map was normalized to the baseline (i.e. T1-T0, T2-T0, and T3-T0), suppression (blue color) in quality and location of perceived sensation is evident in the TENS group compared to no or enhanced sensation in the Sham group.

Functional connectivity
The calculated mean PLI values are presented in figure 4. While Shapiro-Wilk test results revealed that PLI values are not normally distributed, figure 4 shows that the mean PLI value for connections in the TENS group is more than the Sham group in all frequency bands. However, using the non-parametric Friedman test, statistical significance was observed only in two pairs of connections in the gamma-band when comparing TENS and Sham groups. A significant difference was observed in the FC in gamma-band for SI_R-ACC and SI_R-mPFC connections across time phases in the TENS (p = 0.0463 for SI_R-ACC and p = 0.0292 for SI_R-mPFC). However, the same observation was not made in the Sham group (p = 0.956 for SI_R-ACC and p = 0.173 for SI_R-mPFC). In addition to the above, a Mann-Whitney test then examined the lasting effect of these changes (within-group). Results demonstrate that the alteration in SI_R-ACC (T0-T1: p = 0.014, T0-T2: The mean PLI matrix in the gamma frequency band are shown in figure 5(B). This figure shows the functional PLI-based brain network of the gammaband. The size and color of a node represent the normalized sum of the connectivity values related to each node (nodal strength), and the color of pair-wise connections indicates the normalized average PLI value across subjects. Only edges with a normalized mean PLI value >30% (as compared to the overall maximum value) were included in the figure for the sake of simplicity.

Graph analysis
The effect of the TENS intervention on the topology characteristics of the network was examined. Therefore, changes from T0 to T1 in global and local indexes of networks were compared between TENS and Sham groups for the six frequency bands. For both global and local indexes, clustering coefficient, efficacy, and strength were evaluated. Please see figure 6(A). for an illustration of the mean value of extracted global graph measures for the TENS and sham groups. The gamma-band results showed that the global node strength, global clustering coefficient, and global efficiency in the TENS group are higher than those under sham intervention. Results from the independent sample t-test revealed a significant difference between the TENS and Sham groups for the gamma-band network in the global node strength (p = 0.0321), global clustering coefficient (p = 0.0243), and global efficiency (p = 0.0390).
Statistical analysis of the network also revealed a significant difference between the TENS and Sham groups for all local measures in the gamma frequency band (see figure 6(B)) Statistically significant differences in local efficiency were observed, with significantly higher node strength in the contralateral SI (p = 0.0214), contralateral insula (p = 0.0419), ACC (p = 0.0175), and mPFC (p = 0.011) in the TENS group compared with the sham group.
The local cluster coefficients, however, showed no significant effect of the TENS intervention on nodes contralateral to the stimulation. At the same time, the local cluster coefficient difference between groups for ACC (p = 0.0204), mPFC (p = 0.0257), and ipsilateral SI (p = 0.0243) ROIs were found significant. Simultaneously, the independent sample t-test revealed a statistically significant difference in nodal strength indexes for ACC (p = 0.00487) and ipsilateral SI (p = 0.0331) ROIs between TENS and sham groups.

Discussion
The present study showed that functional brain connectivity following TENS intervention significantly differs from those subjected to a sham intervention. Here we investigate the topographical changes in Figure 5. (A) The eight selected source brain regions overlaid on the inflated surface maps using the BrainNet Viewer toolbox (RRID:SCR_009446; www.nitrc.org/projects/bnv/ (B) The normalized, mean FC matrices calculated for the TENS (first row) and Sham (second row) groups at the four different time phases (four columns) in the gamma frequency band. The size and color of a node represent the normalized sum of the connectivity values related to each node (nodal strength) and the color of pair-wise connections indicates the normalized average PLI value across subjects. For the brain network simplicity, the normalized edges with value lower than 0.3 are not depicted. Connections with significant changes (i.e. SI_R to ACC, and SI_R to mPFC) are highlighted with * . Figure 6. Global and local characteristics of functional brain network connectivity. (A) The radar plot depicts the changes (% change in T1 compared to T0) for global indexes as efficiency, cluster coefficient, and strength in six frequency bands on the network extracted by PLI. The red area represents the values for the subjects in the TENS group, while the gray area demonstrates the values in the sham group. While we compared the percent of change, the axes run from the center (−50%) to the outside (+50%) for global measures and run from the center (−10%) to the outside (+10%) for local measures. Independent sample t-test revealed a significant difference in all three global indexes in the gamma band. (B) The radar plot shows the three local index changes for eight nodes in the gamma band. Statistical analysis shows a significant change in SI and Insula contralateral to the stimulation site, ACC, and mPFC in local efficiency between the groups. The local cluster coefficient in ACC, mPFC, and SI ipsilateral to the stimulation site were changed significantly. Moreover, our results revealed that TENS was significantly affected the nodal strength in ACC and SI ipsilateral to the stimulation site. functional brain connectivity and sensory perception alterations associated with TENS intervention.

Effect of TENS on the perception response
The TENS intervention on the evoked perception has shown a significant effect not only between groups but also over time. We have used the VAS scale to rate the perceived sensation following each block of SEPs. Although these measurements are subjective and have limited reliability, we have recorded the location of the subject's evoked sensation as supplementary measurements. We have shown that both cortical and sensation responses are suppressed following TENS. Further experiment with recording sensation following each stimulation impulse is needed to correlate the suppression of cortical activity and perceived sensation responses.

Effect of TENS on functional connectivity
The TENS intervention had a significant impact on gamma-band phase-based FC in the SI-ACC and SI-mPFC pair-wise connections in the contralateral hemisphere to the electrical stimulation. The enhancement of gamma-band FC in the two mentioned connections was remained significant for at least up to an hour following the intervention.
This study illustrated that TENS alters the FC between SI-ACC and SI-mPFC in the gamma band, which occur simultaneously with a decrease in perceived sensation. The aforementioned results can improve our fundamental understanding of the neurophysiological effects on perception and potential association with pain relief mechanisms following TENS intervention. Our results suggest the increased FC between SI-ACC and SI-mPFC in the gamma band may be a possible biomarker for analgesic effect influenced by TENS.
Several pain studies have reported the SI as a termination node for canonical ascending nociceptive pathways, which depict the sensory-discriminative component of pain [45,46]. The gamma band activity of SI has also been reported as an important index for sensory processing [47]. Moreover, ACC has been reported to represent the affective aspect of pain experiences [48]. Additionally, the increase in the ACC's aversive responses and the plasticity in both the SI and ACC in chronic pain patients have been reported previously [49,50]. The anatomical and functional connections between SI and ACC have been previously investigated using optogenetics and fluorescent tracing methods. Singh et al improved the understanding of the sensory processing and neuropathic pain at the cortical level by reporting a direct projection from the SI to the ACC following noxious stimuli [51]. Additionally, gamma oscillations have been proved as an index to represent neural activity in brain FC [17]. Accordingly, our finding on the enhancement of SI-ACC FC in the gammaband following TENS could be possible evidence of the important role of SI-ACC connection for the analgesic effect of TENS.

Effect of TENS on global characteristics of brain network
This study showed an alteration in both global and local brain network metrics in the TENS group compared to the sham group. Using functional MRI, some changes of global metrics of a functional network at infra low frequencies (<0.1 Hz) have been reported in chronic pain patients [52]. Furthermore, in a recent EEG study, Ta Dinh et al demonstrated a gammaband global network reorganization in chronic pain patients [17]. They revealed that chronic pain is associated with a significant decrease in the global efficiency of the brain network. In line with this, our results depict that all three global network parameters (i.e. global node strength, global efficiency, and global cluster coefficient) were generally increased in all frequency bands following TENS intervention. However, results show a significant change only in gammaband.
Moreover, chronic pain has been shown to suppress the global metrics of sensory processing FC networks in gamma-band [17]. As our results showed TENS improve these metrics, the present study suggests that the global nodal strength, global efficiency, and global cluster coefficient could be a potential biomarker to investigate the mechanism of TENS in pain reduction.

Effect of TENS on local characteristics of brain network
In addition to the topological attributes of the brain network, analyzing the network from local aspects is commonly accepted [17,19]. The local modular attributes represent the performance of information exchange with neighboring or directly connected nodes, and information segregation of the network can be assessed by local efficiency and local cluster coefficient [44]. Our results revealed several significant enhancements of local features. Significant increase in efficiency, cluster coefficient, and ACC's strength suggest this ROI as one of the most affected brain areas by TENS intervention. Hautasaari et al has recently shown that that ACC is not involved in early processing of nociceptive stimulation [53], However, in the long-term sensory processing, our results showed TENS could significantly improve the segregation of information in ACC.
Moreover, figure 6(B). shows TENS intervention increased information transmission efficiency for the contralateral SI, ACC, and insula. Alteration in the insular cortex's activity has been reported in sensory (non-painful tactile [54]), acute pain [55], and chronic pain [56] studies. Previous studies also found an anatomical connection of the insula with multiple sensory, motor, limbic, and association areas [57,58]. At the same time, studies show that ACC plays an essential role in sensory-pain processing, and an alteration in ACC activity has been reported in chronic pain patients [51]. Consequently, our result suggests that the local efficiency enhancement of SI, Insula, mPFC, and ACC is another explanation for the mechanism of analgesic effect of TENS.
The changes in the somatosensory cortex activity induced by TENS on the ipsilateral to the stimulation site have been reported along with the contralateral hemisphere [10]. The importance of bilateral brain area activity has recently been reported in sensory and pain networks [59]. Moreover, several studies suggested that applying the TENS on the intact hand for PLP relief on upper limb amputees is more effective than the ipsilateral hand [60,61]. Suppression of interhemisphere FC in the SI following amputation has also been reported in the fMRI study [62]. Consistent with previous research in TENS effect on hemisphere ipsilateral to the stimulation site, a significant change in node strength and cluster coefficient of ipsilateral SI (SI_L) in the gamma band FC network was observed in our study. Thus, local node strength change in the ipsilateral hemisphere to the stimulation site may represent the possible effect of TENS on PLP relief. A possible explanation for this view can be addressed by the gamma-band oscillation role in pain perception, and evidence on the gamma-band oscillation of SI ipsilateral to the stimulated hand might be caused by the same oscillation on the SI to the contralateral stimulated hand.
Our results demonstrated the effect of TENS on the brain FC network involved in sensory/pain processing in healthy subjects. The method used in the present study for inducing the SEPs is based on recently published results in the literature [11] which is different in terms of pulse width and number of pulses when compared with more classic SEP [63][64][65]. In this regard, the utilized method includes double-pulse stimulation with a pulse width of 500 µs to address more charge of electrical stimulation for surface electrical stimulation (compared to intraepidermal) to induce a sensory stimulation. Addtionally, in this study, we examined our objective in healthy subjects to provide an adequately large and homogeneous subject population. Although chronic pain reduction previously has been associated with changes in sensory/pain-related brain FC, the ability of the present TENS intervention to induce such alteration and analgesic effect on pain should be validated in a future clinical trial.

Conclusion
The application of TENS on pain alleviation has been previously proved on chronic pain patients. The underlying brain mechanism following TENS intervention and associated modulation of perceived sensations is not well known. In the present study, we investigate the influence of TENS on sensation and hypothesize that it can be decoded by analyzing the topographical alterations in the FC network of the central nervous system. The results showed a significant increase in the FC between SI_R, ACC, and mPFC (at least for an hour) following TENS intervention. Moreover, results from brain functional network analysis depicted, for the first time, an increase in global and local network characteristics, with a suppression in the intensity and size of the sensation area perceived in response to a sensory stimulation after TENS.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.