Authors: Yuhao Wang, Ryohei Fukuma, Ben Seymour, Huixiang Yang, Haruhiko Kishima, Takufumi Yanagisawa
Categories: Article, Sensory processing, Perception
Source: Communications Biology
Authors: Yuhao Wang, Ryohei Fukuma, Ben Seymour, Huixiang Yang, Haruhiko Kishima, Takufumi Yanagisawa
Insula activity has often been linked to pain perception, making it a potential target for therapeutic neuromodulation strategies such as neurofeedback. However, it is not known whether insula activity is under cognitive control and, if so, whether this activity is consequently causally related to pain. Here, we conducted a double-blind randomized controlled crossover trial to test the modulation of insula activity and pain thresholds using neurofeedback training. Nineteen healthy subjects underwent neurofeedback training for upmodulation and downmodulation of right insula activity using our magnetoencephalography (MEG)-based brain–machine interface. We observed significant differences in insula activity between the upmodulation and downmodulation training sessions. Furthermore, resting-state insula activity significantly decreased following downmodulation training compared to following upmodulation training. Compared with upmodulation training, downmodulation training was also associated with increased pain thresholds, albeit with no significant interaction effect. These findings show that humans can cognitively modulate insula activity as a potential route to develop therapeutic MEG neurofeedback systems for clinical testing. However, the present findings do not provide direct evidence of a causal link between modulation of insula activity and changes in pain thresholds.
Pain perception is determined not only by external stimuli but also by internal states related to perceptual sensitivity. Various neuroimaging studies have demonstrated that pain perception is regulated by a broad cerebral network^1,2^, which includes the S1, S2, insular cortex, anterior cingulate cortex, prefrontal cortex and thalamus^3^.
Among these regions, the insular cortex serves as a key hub in pain perception, contributing to both nociceptive processing and its cognitive modulation^4–10^. Research in animal models has shown that manipulating neural activity in the insula significantly affects pain perception^11,12^. Multiple human studies have shown that the insular cortex is activated following exposure to noxious stimuli^1,7,13–18^. Additionally, insula activity has been reported to be associated with chronic pain states^19^, such as those experienced by patients with fibromyalgia^20,21^ and complex regional pain syndrome^22^. Notably, direct electrical stimulation of the human anterior insula can increase pain thresholds^23^, suggesting that engaging this region can causally influence pain perception. However, despite these associations, it remains unclear whether insular cortex activity can be deliberately controlled and, if so, whether such control has a direct causal effect on pain perception. In other words, a definitive link between insula modulation and changes in pain experience has yet to be established, thus highlighting the need for studies that directly test insula-focused neuromodulation in humans.
One promising strategy to modulate region-specific neural activity is neurofeedback^24^. In a typical neurofeedback system, neural activity is continuously recorded, processed, and displayed in real-time. Participants interact with the feedback, learning to regulate their brain activity toward a target state, thereby facilitating adaptive neuroplasticity^25^. Previous neurofeedback studies have provided preliminary support for targeting the insula in pain modulation. For example, EEG-based neurofeedback training has shown that changes in insular theta-band oscillatory power correlate with alterations in pain perception^4^. Similarly, real-time fMRI studies reported that modulating blood-oxygen-level dependent (BOLD) signals in the insula can lead to changes in pain processing^26,27^. These prior findings suggest a strong relationship between insula activity and pain modulation. Nevertheless, they do not conclusively demonstrate a direct causal link, as they primarily show correlations or group differences rather than providing proof that an individual’s deliberate control of insular activity drives pain changes. Key questions (i) To what extent is insula activity under cognitive control such that it is susceptible to neurofeedback? (ii) Are changes in insula activity associated with changes in pain perception?
To address these questions, we proposed the following (i) the mean cortical excitability, or root mean square (RMS) of the cortical currents, of insula activity can be downmodulated through neurofeedback; and (ii) downmodulating the mean cortical excitability of insula activity through neurofeedback can increase the pain threshold. Additionally, as an exploratory analysis, we examined whether changes in mean cortical excitability were associated with changes in oscillatory power.
To test our hypothesis, we developed a magnetoencephalography (MEG)-based neurofeedback approach that enables participants to modulate their right anterior insula activity in real time. We then conducted a double-blind, randomized controlled crossover trial. Our system was specifically designed to extract the RMS of the insular cortical current—a measure of mean cortical excitability—from the MEG data. By using the RMS of estimated insula currents as the feedback signal, we targeted the overall amplitude of insular activity rather than focusing on a single frequency band. This approach differs from previous insula neurofeedback studies, which primarily targeted oscillatory power (e.g., theta rhythm) or BOLD activity.
Twenty healthy volunteers were recruited and randomly allocated to two groups; 19 of them (11 males, 8 females; mean age 28.8 years; range 22–58 years; all right-handed) completed this trial. A CONSORT diagram of the subjects throughout the trial process is provided in (Fig. 1).Fig. 1Consolidated Standards of Reporting Trials (CONSORT) flow diagram of this double-blind randomized controlled crossover trial.
The subjects engaged in a series of recording and training sessions (Fig. 2a). We chose the right insula region (Fig. 2b) from the automated anatomical labeling (AAL) atlas^28^ as our region of interest (ROI). First, we recorded the MEG signals of the subjects during a 5-min resting period. Subsequently, we measured their pain threshold and proceeded with neurofeedback training. During the training sessions, the subjects were asked to mentally control the size of a black disk displayed on a screen (Fig. 2b). During the upmodulation training, insula activity was directly proportional to the size of the black disk, which indicated that successful training would result in increased insula activity. On the other hand, during downmodulation training, insula activity was inversely proportional to the size of the black disk, implying that successful training would lead to decreased insula activity. Following training, we measured the subjects’ pain thresholds again and recorded brain activity for another 5-min period while the subjects were in the resting state.Fig. 2Neurofeedback training and experimental design.a Each experiment consisted of the following first, MEG resting-state data were recorded for 5 min. Then, the subject’s pain threshold was measured. Neurofeedback training was performed 2 times for 10 min each. Two types of modulation were used to control the feedback picture, one for each experiment. After the pain threshold was measured, the MEG resting-state data were recorded again for 5 min. b Neurofeedback training. The subjects were instructed to intentionally try to enlarge the size of the black disk on the screen.
The total training duration was 1200 s, corresponding to 6000 time points. To compare insula activity between the two types of training, the data were standardized across the entire training process using the magnitude and variance of insula activity during the initial 60 s (baseline) for each subject, with values deviating by more than three standard deviations treated as artifacts and replaced with the median. Time-series data were smoothed using a moving average filter with a window of 100 data points (corresponding to 20 s). This smoothing was performed to better illustrate the trends of neurofeedback-induced modulation of neural activity. Fig. 3a shows the time course of average standardized insula activity, while Fig. 3b illustrates a significant difference in insula activity between the two training conditions. Numerical source data for Fig. 3 are provided in Supplementary Data 1. Specifically, the standardized insula activity during downmodulation training was significantly lower than that during upmodulation training (t = 2.168, p = 0.044 < 0.05; paired t test, n = 19). These results suggested that neurofeedback training successfully modulated insula activity levels.Fig. 3Modulation of insula activity during neurofeedback training.a Temporal dynamics of mean standardized insula activity during upregulation (blue) and downregulation (red) conditions. The solid lines represent the mean activity across subjects, and the shaded areas indicate the standard error of the mean (SEM). All curves in this figure have been smoothed using a 100-point (20-s) moving average to clearly illustrate the neurofeedback modulation trends. The vertical dashed line at 600 s marks a brief pause between two training sessions. b The boxplot compares the mean insula activity during the entire training process for each subject under the upmodulation and downmodulation conditions. Insula activity during upmodulation training was significantly greater than that during downmodulation training (blue for upmodulation training; red for downmodulation training). Two-tailed paired Student’s t test, *p < 0.05. In the boxplots, centerlines represent medians, the box limits represent the lower and upper quartiles, the whiskers represent 1.5 times the interquartile range, and the circles represent outliers.
To further confirm that neurofeedback successfully modulated insula activity levels, we extracted RMS value of estimated right insula currents from the resting-state brain activity before and after training and compared the changes after each training session and differences between the two training types (Fig. 4). Numerical source data for Fig. 4 are provided in Supplementary Data 1. We found that resting-state insula activity significantly decreased after downmodulation training (t = 3.349, p = 0.004 < 0.025; paired Student’s t test, Bonferroni correction; n = 19). The decrease in resting-state insula activity was significantly greater in the downmodulation experiment than in the upmodulation experiment (t = 2.191, p = 0.042 < 0.05; paired Student’s t test, n = 19). This shows that neurofeedback modulation has a persistent effect on insula activity.Fig. 4Alterations in insula activity before and after neurofeedback.a Mean resting-state insula activity before (pre) and after (post) the upmodulation (blue) and downmodulation (red) training are shown. Paired t-tests were conducted to compare pre and post activity levels within each training condition. b Differences in insula activity (post - pre) for the upmodulation (blue) and downmodulation (red) conditions are shown. A paired t-test was conducted to compare these differences between the two conditions. All the data were analyzed by two-tailed paired Student’s t tests with Bonferroni correction. Asterisks (*) denote statistically significant differences, *p < 0.05, **p < 0.01. In the boxplots, dots represent individual data points, centrelines represent medians, box limits represent lower and upper quartiles, whiskers represent 1.5 times the interquartile range, and circles represent outliers.
To ensure that our feedback was indeed specific to insular activity despite the inherent challenges of reconstructing signals from this deep cortical structure, we performed a whole-brain analysis of source activity during both upmodulation and downmodulation training. Specifically, we estimated cortical current distributions across the entire brain for a 300-s resting period, computed the RMS of these currents, and then applied a 10-mm Gaussian smoothing kernel to the cortical surface. We subsequently evaluated the RMS differences (posttraining minus pretraining) between the two feedback conditions by conducting paired t tests on each vertex (downmodulation vs. upmodulation). As shown in (Fig. 5), the RMS difference decreased in the right insular cortex after downmodulation compared with upmodulation. The region with the most significant t value was the right insular cortex, the target of the feedback, although the differences were not statistically significant after correction for multiple comparisons (false discovery rate (FDR), p > 0.05). The result suggests that our magnetoencephalography (MEG) neurofeedback altered the targeted region, the right insula, most effectively.Fig. 5Paired t test results of the comparison of posttraining changes in the whole brain.The color map threshold was set at uncorrected p < 0.05. The bright blue line marks the boundary of the right insula. N = 19, uncorrected.
We also conducted a statistical analysis of pain thresholds pre- and posttraining (Fig. 6). Numerical source data for Fig. 6 are provided in Supplementary Data 1. There was no significant change after upmodulation training (t = −1.816, p = 0.086 > 0.025; n = 19; paired Student’s t test, Bonferroni correction). The pain threshold significantly increased after the downmodulation training (t = −3.103, p = 0.006 < 0.025; paired Student’s t test, Bonferroni correction; n = 19). However, the changes in pain thresholds after the two types of training were not significantly different (t = −0.474, p = 0.641 > 0.05; n = 19; paired Student’s t test).Fig. 6Alterations in pain thresholds before and after neurofeedback training.a Pain thresholds before (pre) and after (post) the upmodulation (blue) and downmodulation (red) training are shown. Paired t-tests were conducted to compare pre and Pain thresholds within each training condition. b Differences in Pain thresholds (post - pre) for the upmodulation (blue) and downmodulation (red) conditions are shown. A paired t-test was conducted to compare these differences between the two conditions. All the data were analyzed by two-tailed paired Student’s t tests with Bonferroni correction. Asterisks (*) denote statistically significant differences, *p < 0.05, **p < 0.01. In the boxplots, dots represent individual data points, centerlines represent medians, box limits represent lower and upper quartiles, whiskers represent 1.5 times the interquartile range, and circles represent outliers.
During the upmodulation training period, there were no significant changes in pain thresholds. The absence of significant changes does not necessarily suggest ineffectiveness. This could stem from the inherent propensity of the training, regardless of its type, to elevate pain thresholds (i.e., a nonspecific effect of time, manifested as habituation in pain thresholds), as detailed in the supplementary information’s repeated-measures ANOVA (Supplementary Table 1 and Supplementary Fig. 1). This inherent effect could counterbalance or mask the specific outcomes associated with upmodulation training.
We conducted a spectral analysis of the resting-state data to compare the changes in right insula spectral power after two types of training (Fig. 7). No significant differences were found. Numerical source data for Fig. 7 are provided in Supplementary Data 1.Fig. 7Comparing the band power between the pretraining and posttraining resting states.We observed no significant difference (n = 19, two-tailed paired Student’s t test, Bonferroni correction) in the spectral power changes (posttraining minus pretraining) after training. In the boxplots, dots represent individual data points, centerlines represent medians, box limits represent lower and upper quartiles, whiskers represent 1.5 times the interquartile range, and circles represent outliers.
Our double-blind randomized controlled trial demonstrates that MEG-based neurofeedback, employing a cognitive modulation paradigm guided by the RMS of estimated insula currents, effectively reduces insula cortical activity. Notably, this modulation occurred not only during active neurofeedback but was also evident in resting-state measurements before and after training, providing robust evidence that the neurofeedback intervention induced sustained changes in cortical excitability of the insula. Although downmodulation of insula activity was accompanied by a significant increase in pain thresholds, the specificity of this effect compared with the upmodulation control condition was not confirmed. These findings suggest that insula activity is amenable to modulation via MEG-based neurofeedback and highlight its potential as a therapeutic target for pain management through repeated neurofeedback training. However, further evidence is needed. Due to the lack of a significant interaction effect between training type and pain threshold changes, these findings must be interpreted cautiously and do not establish a direct causal relationship between insula modulation and altered pain perception. We acknowledge the limitations of our study, particularly the relatively small sample size and the borderline significance of the results, which may affect the reliability and generalizability of our findings. Future studies with larger sample sizes are needed to validate these observations and further investigate the underlying effects.
Most current neurofeedback research employs EEG and fMRI^29^. EEG provides millisecond-level temporal resolution. However, most EEG-based neurofeedback (EEG-NFB) protocols do not modulate specific indicators of neural activity strength. Instead, they focus on regulating spontaneous brain rhythms, which are primarily defined by their oscillatory frequencies and may not directly correspond to the physiological activity of interest^30^. Due to the limited spatial resolution of EEGs, it is challenging to modulate brain activity in specific regions. Unlike EEG-NFB, fMRI-NFB offers millimeter-level spatial resolution and continues to guide participants in successfully modulating brain activity, as indicated by the blood-oxygen-level-dependent (BOLD) signal^31^. However, the lower temporal resolution of fMRI data poses a challenge to the real-time nature of neurofeedback. Researchers have attempted to combine these modalities by using EEG-fMRI integration to balance their strengths but still face difficulties such as electromagnetic interactions between modalities that introduce artifacts affecting the signal-to-noise ratio, thereby reducing data quality compared to acquiring modalities separately and the challenge of effectively integrating the two modalities^32,33^. In our study, we used the RMS of estimated insula currents with a MEG-based brain-computer interface, achieving both high temporal resolution and relatively high spatial resolution and enabling the modulation of the mean cortical excitability of insula activity. Our study provides an alternative direction for the successful modulation of insula activity to improve neurofeedback methods.
The diverse functions of the insula include assessing the emotional aspects of pain, modulating cognitive processes related to pain anticipation and attention, regulating physiological responses (e.g., heart rate and blood pressure), perceiving internal body states, and adjusting pain thresholds on the basis of the sensory input context^34–39^. The successful modulation of insular activity in our research may impact pain perception across these dimensions. Further explorations of how neurofeedback may modulate pain thresholds could be performed in future studies by focusing on specific functions of the insula, thereby increasing our understanding of and ability to manage pain.
In previous studies, brain spectral power adjustments have predominantly been used to modulate insula activity^40^. While some research has indicated potential alterations in pain discrimination^4^, robust double-blind controlled studies demonstrating that changes in insula brain spectral power within specific frequency bands can modify pain perception or thresholds are lacking. In our study, we did not observe significant changes in spectral power within these specific bands that correlated with alterations in insula activity or changes in pain thresholds, suggesting that the overall neural activity in the insula might be more relevant to pain and could be a more appropriate target for neurofeedback.
Several studies have employed brain stimulation techniques such as TMS to modulate brain activity in deep regions like the insula. However, there is ongoing debate regarding whether these methods can effectively reach such deep areas due to limitations in the depth of stimulation. This controversy highlights the challenges inherent in non-invasive brain stimulation techniques. Nasr et al.^41^ discussed the challenges of traditional non-invasive brain stimulation techniques, particularly their limited precision and potential side effects when targeting deep brain regions. They emphasize the need for more advanced approaches that can overcome these limitations. Combination of neurofeedback with closed-loop approaches^42^, especially those that can reliably target deep brain areas such as ultrasound may offer a promising avenue for more effective and precise interventions. Neurofeedback combining with these advanced methods could potentially overcome the limitations of traditional brain stimulation techniques by providing real-time feedback and more accurate targeting.
Although our research was focused primarily on acute pain induced by thermal stimulation, the insula is known to play a role in both the modulation of pain perception and the development of chronic pain^43^. Neurofeedback treatments that modulate insula activities are therefore a plausible clinical option for patients with chronic pain. Two key steps are required to achieve this goal. First, it is necessary to show that insula activity is susceptible to neurofeedback control, and the findings of the current study suggest that this is the case. The second step is to show that training can have an impact on pain perception. For most neuromodulatory strategies, including neurofeedback, regular, repeated training sessions are required. The increased pain threshold we observed after a single session of neurofeedback is in keeping with this prediction, but for robust demonstration against the control condition, multiple sessions are likely required.
In conclusion, we demonstrated that MEG-based neurofeedback can be used to modulate insula activity, with a potential impact on pain perception. Our approach using estimated cortical currents in the insular cortex offers a distinct perspective on how to modulate brain activity through neurofeedback. Future studies should investigate the long-term effects of repeated neurofeedback training on pain perception and evaluate the efficacy of this approach in clinical populations with chronic pain.
We recruited subjects through online solicitation. A total of 20 healthy volunteers, including 11 males and 9 females, were initially recruited for the study. One female subject withdrew partway through the study, resulting in 19 subjects (11 males and 8 females) completing the trial. We collected information solely on the biological sex of the subjects and did not gather data regarding their cultural gender. Efforts were made to ensure a balanced representation of sex among the participants. The mean age of the subjects was 28.8 years, with an age range of 22–58 years. All the subjects were right-handed, as determined by self-reports.
The study adhered to the Declaration of Helsinki, was performed per protocols approved by the Ethics Committee of Osaka University Clinical Trial Center (No. 18472-2) and was registered in the UMIN Clinical Trials Registry (UMIN000037063). We obtained written informed consent from all subjects after informing them of the purpose and possible consequences of this study. All ethical regulations relevant to human research participants were followed. Subjects received a compensation of 5000 JPY for each experimental session they attended, which was independent of their experimental performance.
Before starting the study, we prepared a randomization list using block randomization with block sizes of 2 and 4. The list was kept confidential by the experiment operator. This list equally divided the twenty subjects we intended to recruit into two groups. The blocked randomization list was incorporated into the experiment initiation program, which was subsequently encrypted. During the setup, the operator input the subject’s ID to launch the neurofeedback program. Through this process, the subjects and the experiment operator were blinded to the training allocation. The data preprocessing and analysis procedures were conducted using participants’ individual IDs, without incorporating any group assignment information.
To assess the pain threshold, contact thermal stimuli were applied using the Medoc Pathways Pain & Sensory Evaluation System (Medoc Ltd., Ramat Yishai, Israel). A thermode was gently applied to the ventral side of the first web space on the subject’s left hand. We employed the limit method, wherein the thermode temperature was elevated at a steady rate of 1.5 °C/s from a starting temperature of 32 °C. The subjects were instructed to press a button when they perceived pain. This measurement was conducted 15 times. The initial 3 attempts were treated as trials and disregarded. The subject’s pain threshold was determined as the average temperature when the subject pressed the button.
MEG signals were measured using a 160-channel whole-head MEG system equipped with coaxial gradiometers (MEGvision NEO, Yokogawa Electric Corporation, Kanazawa, Japan) housed in a magnetically shielded room. A projection screen was used to present visual stimuli (Presentation, Neurobehavioral Systems, Albany, CA) from a liquid crystal projector (LVP-HC6800, Mitsubishi Electric, Tokyo, Japan). The data were captured in a frequency band ranging from 0.1 to 500 Hz and sampled at a rate of 2000 Hz.
Every subject assumed a supine posture within a magnetically protected chamber. Prior to initiating the MEG measurement, 5 fiducial coils were affixed to the subject’s face. The coil positions were ascertained both prior to and following each MEG session to establish the relative orientation and location of the MEG detectors to the subject’s cranial structure. The permissible maximum deviation in coil position from the beginning to the end of the recording was set to 5 mm. To coregister the MEG data with individual magnetic resonance imaging (MRI) scans (T1-weighted; Signa HDxt Excite 3.0 T, GE Healthcare UK Ltd., Buckinghamshire, UK), we captured the three-dimensional facial topology and mapped 50 points across each subject’s scalp using FastSCAN Cobra (Polhemus, Colchester, Vermont, USA). The 3D facial surface data were superimposed on the anatomical facial surface data extracted from the MR images.
During the resting-state recording phase, subjects were instructed to remain awake throughout the MEG scan, keep their gaze fixed on the central fixation point on the screen, and abstain from engaging in any taxing cognitive processes. Under these conditions, continuous MEG data were collected for a span surpassing 300 s for each resting-state session.
During the neurofeedback training phase, the subjects were directed to maintain their focus on the central fixation point on the screen. The subjects were tasked with exerting effort to enlarge the black disk displayed in the center of the screen. The total duration of the neurofeedback training was 1200 s, split into two sessions of 600 s each. Between sessions, we inquired about the subjects’ feelings and strategies. No additional breaks were provided during the training, except for this interval.
For online estimation of the cortical currents in right insular cortex from MEG data, we used Variational Bayesian Multimodal EncephaloGraphy (VBMEG)^44,45^ (ATR Neural Information Analysis Laboratories, Kyoto, Japan) in MATLAB (version R2013). First, cortical segmentation was performed on T1-weighted MR images via FreeSurfer software (Martinos Center for Biomedical Imaging, Massachusetts General Hospital, Charlestown, MA, USA), generating a subject-specific cortical surface model. On this cortical surface, 4004 current dipoles were predefined, each positioned at a fixed location covering the entire cortex and oriented perpendicular to the local cortical surface. VBMEG was then used to estimate the time-varying current amplitude at each dipole by applying an MEG/EEG inverse solution incorporating fMRI-informed priors. To localize the right insular cortex, as defined in the AAL atlas, in each participant’s anatomical MRI, the individual MRI was spatially normalized to the MNI152 template space using SPM8 (Wellcome Trust Centre for Neuroimaging, University College London, UK). This process generated a deformation field and an affine transformation matrix that mapped the subject space to MNI space. Following normalization, the inverse transformation was applied to the AAL atlas, which was originally defined in MNI space, to warp it into each participant’s native anatomical space. The right insular cortex label was then projected onto the subject’s cortical surface model, ensuring anatomical alignment with the individual’s brain structure. The transformed region was used as the region of interest for source localization and neurofeedback training within the VBMEG system. The lead field was computed by applying the boundary element method to the coregistered sensor positions and polygon model based on data collected immediately before each 10-min training session. Cortical currents were estimated from the MEG signals recorded from 110 selected sensors (Supplementary Fig. 3.). The hyperparameters m0 and γ0 were set to 100 and 10, respectively. This online cortical current estimation method was chosen because our previous researchs^25,46^ have confirmed its effectiveness in evaluating motor cortical responses for neurofeedback applications.
Before the neurofeedback training session, MEG signals were first recorded for 300 s without the subject in the scanner. This environmental noise data was used to characterize baseline noise for signal preprocessing and to refine the VBMEG reconstruction filters. The neurofeedback training session began with a 300-s resting-state MEG recording, during which participants remained calm and fixated on the center of a disk presented in front of them. VBMEG was then applied to the resting-state MEG data to estimate the cortical currents at 4004 vertices over the 300-s period (Fig. 2b). Among these, we utilized the cortical currents \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {{{{\rm{I}}}}}_{{{{\rm{v}}}}}\left({{{\rm{t}}}}\right)
\usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \left(v=1\ldots 225\right) $$\end{document}v=1…225 of the right insula cortex, as defined by the AAL2 parcellation provided in VBMEG. We define \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ T $$\end{document}T as the total number of time points in the 300-second baseline period. For the estimated cortical currents of 300-second duration, we calculated the time average of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {{{{\rm{I}}}}}_{{{{\rm{v}}}}}\left({{{\rm{t}}}}\right) $$\end{document}Ivt (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\mu }_{{{{\rm{v}}}}} $$\end{document}μv) and standard deviation (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\sigma }_{{{{\rm{v}}}}} $$\end{document}σv) for each vertex as 1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\mu }_{v}=\frac{1}{T}{\sum }_{t=1}^{T}{I}_{v}\left(t\right),{\sigma }_{v}=\sqrt{\frac{1}{T}{\sum }_{t=1}^{T}{\left({I}_{v}\left(t\right)-{\mu }_{v}\right)}^{2}} $$\end{document}μv=1T∑t=1TIvt,σv=1T∑t=1TIvt−μv2 Then, using \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {{{{\rm{\mu }}}}}_{{{{\rm{v}}}}} $$\end{document}μv and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {{{{\rm{\sigma }}}}}_{{{{\rm{v}}}}} $$\end{document}σv, we obtained the standardized current \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \hat{{{{{\rm{I}}}}}_{{{{\rm{v}}}}}}\left({{{\rm{t}}}}\right) $$\end{document}Iv^t for each v-th vertex at time t as 2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \hat{{I}_{v}}\left(t\right)=\frac{{I}_{v}\left(t\right)-{\mu }_{v}}{{\sigma }_{v}} $$\end{document}Iv^t=Ivt−μvσv The standardized current \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \hat{{I}_{v}}\left(t\right) $$\end{document}Iv^t of 300 s was segmented into 300 intervals of 1000 ms without overlaps, within which the standardized currents were time averaged to compute the mean values for the *k*-th interval, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \bar{{I}_{v}}\left({{{\rm{k}}}}\right) $$\end{document}Iv¯k. We define \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {{{\rm{V}}}} $$\end{document}V as the number of vertices in the right insular cortex. Then, we computed the RMS value of all vertices within the right insula cortex for the k-th time interval as 3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ A\left(k\right)=\sqrt{\frac{1}{V}{\sum }_{v=1}^{V}{\left(\bar{{I}_{v}}\left(k\right)\right)}^{2}} $$\end{document}Ak=1V∑v=1VIv¯k2 The baseline insula activity and its standard deviation across 300 intervals were calculated as 4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\mu }_{A}=\frac{1}{300}{\sum }_{k=1}^{300}A\left(k\right),{\sigma }_{A}=\sqrt{\frac{1}{300}{\sum }_{k=1}^{300}{\left(A\left(k\right)-{\mu }_{A}\right)}^{2}} $$\end{document}μA=1300∑k=1300Ak,σA=1300∑k=1300Ak−μA2 Similarly, we averaged the real-time estimated currents of the *v*-th vertex in the right insular cortex using a 200-ms time window and subsequently standardized the data by using \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\mu }_{v} $$\end{document}μv and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\sigma }_{v} $$\end{document}σv to compute the real-time insula activity. The real-time RMS value for insula activity was calculated as 5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \hat{A}\left(t\right)=\sqrt{\frac{1}{V}{\sum }_{v=1}^{V}{\left(\hat{{I}_{v}}\left(t\right)\right)}^{2}} $$\end{document}A^t=1V∑v=1VIv^t2 During neurofeedback training, the feedback signal was computed differently for upmodulation and downmodulation sessions on the basis of real-time insula activity. Upmodulation Feedback: For upmodulation, the aim of feedback was to increase insula activity. The feedback value was calculated as 6\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {{{{{\rm{out}}}}}_{-}{\mbox{vals}}}_{{\mbox{up}}}=\left(\frac{\hat{A}\left(t\right)-{\mu }_{A}}{8\cdot {\sigma }_{A}}+0.5\right)\times 255 $$\end{document}out−valsup=A^t−μA8⋅σA+0.5×255 Downmodulation Feedback: For downmodulation, the aim of feedback was to decrease insula activity. The feedback value was calculated as 7\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {{{{{\rm{out}}}}}_{-}{\mbox{vals}}}_{{\mbox{down}}}=\left(-\frac{\hat{A}\left(t\right)-{\mu }_{A}}{8\cdot {\sigma }_{A}}+0.5\right)\times 255 $$\end{document}out−valsdown=−A^t−μA8⋅σA+0.5×255 Clipping and Rounding: The feedback values out_vals were rounded to the nearest integer and clipped to the range [0, 255]:8\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {{{{\rm{out}}}}}_{-}{\mbox{vals}}=\min \left(255,\max \left(0,{\mbox{round}}\left({{{{\rm{out}}}}}_{-}{\mbox{vals}}\right)\right)\right) $$\end{document}out−vals=min255,max0,roundout−vals Notably, we did not remove artifacts to select the noisy time segments for the calculation for neurofeedback training. All the MEG signals obtained in the resting-state and neurofeedback sessions were used for the calculations above. For both the upmodulation and the downmodulation conditions, the subjects were instructed to intentionally try to increase the size of the black disk on the screen. The neurofeedback strategies employed by participants and their perceived success rates in modulation were documented at the end of the experiment through unstructured interviews. We selected the RMS of the insular source signals as the neurofeedback metric (instead of, for example, a frequency-domain power measure spanning 1–100 Hz) because the RMS provides a direct and efficient index of overall insular cortical activity in real time. Mathematically, the RMS of a broadband signal is proportional to the total power of that signal across the defined frequency range, meaning that the time-domain RMS effectively captures the aggregate neural power without requiring an explicit spectral decomposition. This approach takes advantage of the high temporal resolution of MEG by allowing immediate computation from the most recent data, without relying on sliding time windows or frequency transforms, thereby ensuring minimal latency and computational overhead. In practice, the RMS-based feedback is also straightforward to an increase in the RMS value corresponds directly to an increase in the insula’s signal amplitude (and thus overall neural activation), making it an intuitive measure of cortical excitability for both experimenters and participants. However, we acknowledge that using the RMS as the feedback signal lacks frequency specificity. Because the RMS collapses activity across all frequencies, it cannot distinguish which particular oscillatory bands (e.g., theta, alpha, beta, gamma) contribute to changes in insula activity. In other words, any frequency-specific neural dynamics within the insula are obscured when using a broadband amplitude metric. This trade-off was deemed acceptable for our purposes, given that our goal was to modulate the general level of insular activity rather than target a specific EEG or MEG frequency band. We opted not to use a frequency-domain power calculation for the feedback (such as summing the power spectral density from 1 to 100 Hz in each interval) because it offered no clear advantage over the time-domain RMS while incurring additional complexity. A frequency-domain approach would require computing power spectra for each time segment and integrating across frequencies, operations that could introduce extra latency or noise in the real-time feedback loop. In contrast, the RMS measure yields an equivalent reflection of total insular activity power in a simpler and more time-responsive manner, albeit at the expense of detailed spectral information. ### Resting-state data processing Resting-state MEG data processing was performed using Brainstorm, which is an open-source application available under the GNU general public license ([http://neuroimage.usc.edu/brainstorm](http://neuroimage.usc.edu/brainstorm)), and MATLAB version R2017b. To ensure consistency with the online analysis, we also utilized the selected 110 channels. Cortical surfaces reconstructed from each subject’s T1-weighted anatomical MRI scan using FreeSurfer were imported into Brainstorm. The standard preprocessing pipeline was applied; the data were first visually inspected to remove artifacts, such as eye blinks and muscle movements. Artifacts were removed through visual inspection and by applying the Signal Space Separation method. Cardiac and ocular artifacts were further minimized by detecting and removing signal segments corresponding to heartbeats and eye movements using event-based artifact rejection techniques. Bad channel detection and replacement, and head position correction; then, we applied a notch filter at powerline frequencies and harmonics to attenuate powerline artifacts. A bandpass filter was used to remove low-frequency signals (<0.3 Hz). Source reconstruction was performed using the minimum-norm estimation (MNE) method, We applied the dynamic statistical parametric mapping (dSPM)^47^, which is derived from the default MNE and standardized via the noise covariance matrix. The resulting dSPM maps represent a set of z scores, expressed in unitless “z” values. Therefore, we used the noise covariance matrix, which was generated from the MEG signals of 5 min of empty-room data, to the source reconstruction in Brainstorm. Then, the individual source was projected into default MNI ICBM152/FreeSurfer 15000 vertices anatomy. The time course of the source activity was extracted from the right insular cortex using the AAL2 template provided in Brainstorm. The RMS value of the mean absolute estimated currents across all vertices within the right insula cortex was extracted to represent resting-state insula activity. The relative power spectral density (PSD) was estimated from source activity using the Welch method. We computed the average PSD power in the delta (0.3–4), theta (4–8 Hz), alpha (8–13 Hz), beta (13–30 Hz), and gamma (30–100 Hz) bands. It should be noted that we used two source reconstruction methods to ensure reproducibility and achieve high-performance neurofeedback. MNE in Brainstorm was applied for offline analysis of MEG signals during the resting state, while VBMEG was used for online neurofeedback experiments. MNE was chosen for offline analysis because it is a commonly used method, facilitating easier reproduction of our results compared to VBMEG. On the other hand, VBMEG can estimate cortical currents with high accuracy in real time, and its efficacy has been validated in our previous studies^25,46^. therefore, we constructed our neurofeedback system with VBMEG. ### Insular activity The insular activity is quantified by computing the RMS value of estimated currents across all vertices within the right insula cortex for each time segment. This RMS value represents the level of insula activity during that specific interval. The formulas are provided below, with symbols defined 9\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ A\left(k\right)=\sqrt{\frac{1}{V}{\sum }_{v=1}^{V}{\left(\bar{{I}_{v}}\left(k\right)\right)}^{2}} $$\end{document}Ak=1V∑v=1VIv¯k2 ### Sample size calculation We used G*Power 3.1 to calculate the required sample size for a paired *t* test for within-subject comparisons, specifically for changes in insula activity during training and pre/post changes in insula activity and pain thresholds. The effect size (Cohen’s d^48^) was set to 0.7 (large effect), with a power of 0.8 and a significance level of 0.05. The results indicated that at least 19 participants were required for the study. Based on this calculation, we initially set the sample size to 20 participants. Although relatively small, this sample size (*n* = 20) is consistent with previous fMRI and MEG decoded neurofeedback studies, which typically include 10–20 participants^25,27,46,49–54^. Post hoc effect size calculations indicate that insula activity measures during neurofeedback training and resting-state changes exhibited moderate effect sizes (0.497 and 0.503, respectively). However, the effect size for pain threshold changes was substantially smaller (0.109), suggesting that our initial estimate of 0.7 may have been overly optimistic. To improve the reliability of our findings, we plan to increase the sample size in future studies. ### Statistical analysis and reproducibility Statistical analysis was performed via Python 3.8.5 with the scipy package v.1.5.2^55^. Two-sided statistical tests were performed for all the subjects. To test Hypothesis 1 (pre- versus posttraining changes in insula activity), we performed a paired *t* test to compare insula activity measured before and after neurofeedback training. To evaluate Hypothesis 2 (changes in pain thresholds associated with modulation of insula activity), we similarly employed a paired *t* test on pain threshold data measured before and after training. Exploratory analyses assessing oscillatory power changes across multiple frequency bands were also conducted using paired *t* tests. To control for multiple comparisons arising from related statistical tests (e.g., analyses of multiple frequency bands in exploratory testing), we applied Bonferroni correction by adjusting the original significance threshold (*α* = 0.05) according to the number of comparisons performed (n). Specifically, the adjusted threshold for statistical significance was set as α = 0.05/n for each set of related comparisons. ### Reporting Summary Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article. ## Supplementary information Transparent Peer Review file Supplementary information Description of Additional Supplementary Files Supplementary data 1 NR-reporting summary