With the development of functional neuroimaging, our understanding of pain has matured to a concept of multidimensional experience in which the brain integrates inputs from sensory-discriminative and cognitive-emotional systems as a central hub [1]. Pain neuroimaging has also proved that complementary medicine approaches, beyond pharmacological analgesic means, can modulate these central systems [2].

Mindful breathing (MB) is widely accepted as an authentic treatment for pain relief by patients and society in general [3]. The adoption of MB is a welcomed change in our clinical mindset. It decreases our tendency to rely exclusively on pain medications, which can sometimes escalate to dire side effects [4]. In addition, MB techniques are self-facilitated and easy to implement compared with other methods. MB requires learners to regulate their attention to the dynamic interoceptive nature of breathing. When other thoughts disrupt the focus, the learners need to recognize the disruption and refocus on their breathing. This practice gradually gains a learner’s mental control and stabilization abilities. This ability has been proven to alleviate anxiety, stress, depression, and pain among patients [5].

However, the brain mechanisms for the MB practice in pain modulation are poorly understood [5]. Pain is believed to be represented in the brain via affective and sensory networks [6]. Briefly, the ascending noxious signal reaches the spinal trigeminal nucleus, thalamus, and sensory cortex [7,8]. The signal is also processed in the insular cortex and subjectively evaluated in the anterior cingulate cortex, prefrontal cortex (PFC), and other cognitive-emotional regions [9,10]. However, whether the two dimensions can be separately modulated is not well supported by the existing literature [2]. One study investigated a group of long-term Zen meditation practitioners with significantly higher pain thresholds and found increased activation in sensory-related regions (thalamus and insula) but reduced activation in pain-evaluation areas (medial PFC and anterior PFC [aPFC]) [11]. Moreover, direct associations were found between the level of PFC deactivation and meditation-induced pain reduction [11]. In addition, another study compared the effect of a real-meditation training program and a placebo relaxation program, which found increased connectivity between the posterior cingulate cortex and the dorsolateral PFC (DLPFC) [12,13]. Collectively, these findings suggest that mindfulness meditation may modulate pain through a unique mechanism (eg, high-level cortical function in the PFC).

An existing problem with current MB training is the difficulty of long-time attention focusing, especially for beginners, as this is a subjective interoception process. Therefore, there is an urgent need for a tangible method that can provide an immersive sensory guide to facilitate mental control and match the expectations of the current tech-savvy generation. Recently, the development of virtual reality (VR) has enabled the implementation of such methods. VR is a computer-simulated 3D and interactive experience [14]. Delivered by a visual-audio headset, the virtual experience modulates human sensory and emotional systems. The VR technology has effectively managed pain from burns, cancer, and dental procedures [15,16]. Researchers assumed that this process translocates the patients into immersive 3D auditory (eg, esoteric music and rain sound) and visual (eg, breaking waves and moving geometric patterns) contextual experiences [17], although the exact brain mechanisms remain unclear. For a long time, the VR-based pain modulation effect was understood as a distraction mechanism, in which limited attention is partially occupied by exteroceptive VR stimulation (auditory and visual) instead of pain [15,16]. A couple of previous studies found decreased activation of pain-related regions during VR sessions. However, recent studies have shown an after-VR effect rather than only distraction mechanisms [18].

Although proven to be effective in pain modulation, the underlying brain mechanisms of the two processes—abstract sensory-interoception and VR-based sensory-exteroception—remain unclear. This study compared the effectiveness of the two methods in modulating the patients’ pain thresholds in the same study design and clinical environment. We used a week-long protocol in which the participants practiced the traditional MB (TMB) and VR breathing (VRB) in the lab on the first and seventh day, respectively, intercalated by five daily MB practices at home. During the in-lab sessions, we measured the participants’ pain threshold using a facial thermal quantitative sensory test (tQST) after their breathing practices. We used functional near-infrared spectroscopy (fNIRS) as a neuroimaging technique to measure participants’ cortical connectivity and activation. fNIRS is a novel optical brain imaging technique that uses near-infrared light to monitor oxygen levels at multiple cortical locations [19]. Although with less spatial resolution and limited light penetration ability, studies have found that fNIRS signals are highly correlated with blood-oxygen-level-dependent signals [20]. It can be a promising substitute for fMRI in many particular scenarios, given that it is quiet, nonferromagnetic, and relatively motion tolerant. In addition, the higher resolution of fNIRS can provide more physiological information, such as heart rate variability (HRV) [21,22]. In the field of pain, an increasing number of fNIRS-based investigations have been developed [23,24]. In this study, owing to the noise level requirements of mindfulness meditation, and the use of electronic VR devices, we selected fNIRS as the neuroimaging technique to study the VRB and TMB practices by focusing on the key cortical regions for sensory processing and inhibition.

We recruited 40 healthy adult participants in this study (women: 21/40, 52%; age: mean 28 years, SD 4 years). Our exclusion criteria included significant hearing or visual impairment, a history of chronic pain or recent acute pain, significant medical conditions, or current evidence of respiratory distress or asthma. The recruited participants were randomly divided into 2 groups. The first group had interoceptive breathing focusing sessions using an in-house–developed, visual-auditory, 3D VR technology aid (VRB: n=20; women: 11/20, 55%; age: mean 26 years, SD 4 years), whereas the second group had abstract MB (TMB: n=20; women: 10/20, 50%; age: mean 29 years, SD 4 years). This study was approved by the institutional review board of the University of Michigan.

We designed a week-long protocol for both groups. Within the seven days of the protocol, we scheduled each participant for two in-person appointments in the lab on the first and seventh days. From the second to the sixth day, we asked the participants to practice home self-guided exercises following the instructions (Multimedia Appendix 1).

For both in-lab sessions, we asked each participant to complete a set of McGill Pain Questionnaire and the Positive and Negative Affect Schedule (PANAS) questionnaire. Participants were then seated in a dental chair, and we helped them with the comfortable placement of headphones, Oculus Rift virtual imaging equipment (Oculus VR), and a plethysmography belt. Finally, we set up the fNIRS imaging sensors and thermal quantitative sensory thermode, as indicated in Figure 1 (A and B).

As indicated in Figure 1C, we first asked the participants to relax and rest in the dental chair for 5 minutes while we collected the resting-state fNIRS data. The participants then underwent a 10-minute interoceptive breathing awareness practice (TMB or VRB protocol). Participants in the VRB group watched an in-house developed VR display of a 3D lungs image that was synchronized to their inhaling and exhaling cycles in real time via an Oculus Rift device (Oculus VR), as shown in Figure 1B. Meanwhile, they listened to their breathing sounds using headphones. In an early study, Abushakra and Faezipour [25] conceptualized a mobile app with synchronization based on the breath sounds picked by the microphone; however, in practice, their concept could not accurately discern between expiration (breathing out) and inspiration (breathing in) sounds. In our independent method, breathing synchronization was performed using a Braebon plethysmography belt (Great Lakes Neuro Technology). Alternatively, we asked the TMB group participants to abstractly imagine their breathing inflating and deflating. We collected fNIRS data during breathing practice for both groups.

Next, we administered 20 trials of the tQST (Medoc Pathway System). We used 20 times repeated measurements and used the averaged temperature thresholds for further analysis to ensure the test-retest. We placed a single unilateral thermode on the left mandibular nerve branch of the trigeminal cranial nerve (V3 division) for each participant. Within each trial of the 20 trials, the thermode temperature controlled by the controlling device increased from a baseline of 30 °C (86 °F) to a maximum temperature of 50 °C (122 °F), with an increase rate of 1 °C per second. We instructed participants to click the button on the mouse at the first detection of pain, as it stopped the temperature from increasing. The thermode temperature then returned to its baseline and gave the subject a 10-second rest period before the next thermal trial. We also collected fNIRS brain data during the tQST session.

Upon completion of all 20 tQST trials, we asked participants to relax for another 5 minutes for a final collection of resting-state fNIRS data. In addition, we asked participants to complete another set of PANAS and McGill pain questionnaires. Following the completion of the first session, we gave participants a sheet of at-home breathing practice prompt and instructed them to read and complete this exercise for 5 minutes three times a day (after waking, midday, and before bed). Finally, we asked participants to repeat the same protocol during their second in-lab visit on day 7.

We used a TechEN-CW6 fNIRS system (TechEn, Inc) with wavelengths of 690 and 830 nm. The fNIRS cap setup included eight emitters of near-infrared light and 28 detectors spaced 3 cm apart, yielding 45 data channels (CHs) deployed at the bilateral aPFC, premotor cortex (PMC), supplementary motor area (SMA), motor cortex, primary somatosensory cortex (S1), and visual cortex (V1), as indicated in Figure 2. Neuroimaging data were collected at a sampling rate of 25 Hz throughout the entire experiment.

The probe holding cap was established and applied consistently for each participant using the international 10-10 transcranial system positioning [26]. We designed the cap in three sizes—56, 58, and 60, respectively (an example cap is shown in Figure 2), to account for head size variation. In addition, we applied a photogrammetry method to register all optodes and data CHs onto the cortical surface. The detailed method was described in our previous paper [27]. Briefly, we used the Structure Sensor (Occipital Inc) with an iPad (Apple Inc) to capture the 3D photos of the designed caps in three sizes. We then loaded the 3D photo in the MATLAB software (Mathworks) and pinpointed the locations of fNIRS optodes with five fiducial markers (Nasion, Inion, Cz, AR, and AL in the 10-10 system). The derived optodes coordinates were affinely transferred into the Montreal Neurological Institute space using the MATLAB-based AtlasViewerGUI toolbox (Citation). The midpoints between the source and detector (optodes) pairs were used as the coordinates for each CH. Finally, we matched the regions detected by each CH using the estimated center points in the neurosynth.org database. We also estimated the covering range for each CH with a voxel size of 10 mm using the WFU_pick atlas in the XJview toolbox [28].

We estimated the HRV data using time-domain methods by calculating the SD of the normal-to-normal parameter using MATLAB software (MathWorks). We calculated the parameters based on the optical density data of 830 wavelengths, bandpass filtered with a cutoff frequency of 0.01-2 Hz. The parameters were then calculated using the formulas also available in the paper by Wang and Huang [29]:

where N is the total number of R peaks and is the mean of the R-R intervals.

We analyzed fNIRS data using the near-infrared spectroscopy–toolbox [30] in MATLAB software (MathWorks) and a set of customized scripts. In this study, we focused only on oxyhemoglobin (HbO) and not deoxyhemoglobin (HbR). Quantitative analysis indicated that HbR signal changes contributed 16%-22%, whereas HbO signal changes contributed 73%-79% to the total changes measured by fNIRS in cortical hemoglobin concentrations [31].

We applied a generalized linear model with prewhitening and robust least squares [32] to analyze the data collected during the tQST session. Specifically, the raw fNIRS data were first down-sampled to 2 Hz and then converted into HbO and HbR using the modified Beer-Lambert law [33]. We then applied a CH-based generalized linear model regression to each participant’s data, assuming a canonical hemodynamic response function model peaking at 6 seconds. The process can be expressed as follows:

where W is the whitening matrix, y is the observed HbO data, x is the design matrix (model), β is the regression coefficient, ε is the residual, and i and j separately represent the participant and CH index, respectively.

Group-level analysis was conducted using a linear mixed-effects model based on the regression coefficients derived from the individual-level analysis. The model can be expressed as follows:

where Y_g is the regression coefficient obtained from the first level, X_g is the fixed effects term including the modeled brain response at the group level, Z_g is the random effects term counting for between-participant difference, B and θ are the fixed and random effects (coefficients) at the group level, and ε is the residual. Finally, we used a 2-tailed t test to examine the effect of a specific CH:

where Cov_group is the covariance of the group-level model.

To study the brain mechanism during breathing practice, we calculated the associated functional connectivity patterns using the pipelines in the near-infrared spectroscopy–toolbox [30]. The calculation process was described in a previous study [34]. The raw fNIRS data were first down-sampled to 4 Hz. We converted the raw data into HbO and HbR using the modified Beer-Lambert law [33]. We then used bandpass filters to filter the HbO data into two frequency bands: high (0.5-1 Hz) and low (0.01-0.08 Hz) frequency bands. These two frequency bands were selected to avoid the physiological signal bandwidth, including the Mayer wave (0.1 Hz), respiratory (0.3-0.5 Hz), and cardiac (1-1.5 Hz) relevant fluctuations [35]. Next, we calculated the between-CH correlation at the individual level using the robust correlation method in the toolbox [34]. Then, the individual-level correlation coefficient was converted to a Z score using Fisher Z-transform [36]. Finally, a linear mixed-effect model was applied to obtain the group-level connectivity effect. This calculation was implemented on the data collected from both groups during the two lab visits.

We inspected the relationship between brain connectivity during breathing practice and the temperature threshold measured during the tQST session. We ran an elastic net regression to select the best region-to-region connections for the temperature thresholds. This selection process was performed using the Lasso toolbox in MATLAB (MathWorks). Specifically, we iteratively varied the weight controlling the lasso versus ridge optimization from 0 to 1 in increments of 0.1 to achieve a minimum squared error. With each weight, we applied 10-fold cross-validation with Monte Carlo repetitions (100 times) to guarantee a converged output. Next, we calculated the Pearson correlation coefficient (r) and Spearman rank correlation coefficient (ρ) between the selected connections and temperature thresholds. In addition to the Pearson correlation coefficient, we calculated the Spearman rank correlation coefficients to prevent outliers from affecting the correlation analysis [37].

Our first observation was that the average temperature threshold measured by the tQST increased from 45.4 to 46.0 °C (P=.001) for the TMB group, whereas the threshold increased from 46.5 to 47.1 for the VRB group (P=.02) from in-lab sessions 1 to 2 after one week (Figure 3A). Although the thresholds increased at approximately the same level, we did not find a between-group difference. In addition, Figure 3 (B and C) indicate that participants in the TMB group gained higher serenity scores after the practice (visit 1: P=.01; visit 2: P=.001). However, they felt more tired (visit 1: P=.03; visit 2: P=.01), according to the PANAS questionnaires.

Figure 4 shows the HRV estimation results, where panel A shows the heartbeat signal extracted from the fNIRS signal, and panel B shows the estimated SD of normal-to-normal parameters for different sessions. We did not find significant differences between the TMB and VRB groups, suggesting that there might be no different breathing patterns associated with the two breathing practices.

To study and compare the mechanisms of the two breathing practices, we then investigated the brain connectivity during both breathing practices (pre-tQST) and brain activation during the tQST sessions (thermal pain challenge). As shown in Figure 5A, during breathing practices, the TMB group demonstrated a denser functional connectivity pattern (P<.001) among areas such as the aPFC, PMC, SMA, the S1, and the auditory and visual regions than the VRB group.

We further investigated the correlation between brain connectivity and temperature thresholds across participants. As shown in Figure 5 (B and C), the connections within the aPFC (r=−0.46, P=.003; ρ=−0.52, P<.001), between the bilateral PMC (r=−0.49, P=.001; ρ=−0.45, P=.004), between the aPFC and PMC (r=−0.47, P=.002; ρ=0.32, P=.046), and between the S1 and PMC/SMA (r=+0.47,0.48, P=.002; ρ=0.52, 0.51, P<.001) were found to be associated with the temperature thresholds in the TMB group. Whereas in the VRB group, the connections between the PMC/DLPFC and the inferior parietal lobe (r=−0.51, P<.001; ρ=−0.50, P=.002), between the S1 and V1 (r=−0.51, P=.002; ρ=−0.42, P=.01), and between the S1 and superior temporal gyrus (STG; r=−0.48, P=.003; ρ=−0.41, P=.01) were found to be associated with the temperature thresholds in the VRB group.

Next, we examined cortical activation during the tQST session to study how the brain processes pain after the two types of practices. As indicated in Figure 6, the analysis first confirmed contralateral S1 region activation in both groups at both visits (TMB group visit 1: CH 35, t₇₂=3.0, P=.004; TMB group visit 2: CH 35, t₇₂=2.1, P=.04; VRB group visit 1: CH 35, t₇₂=2.4, P=.02; CH 30, t₇₂=2.4 P=.02; CH 29, t₇₂=4.9, P<.001; VRB group visit 2: CH 30, t₇₂=3.4 P<.001).

Our second observation was the activation of the aPFC region in both groups, as indicated in Figure 6 (TMB group, visit 1: CH 2-8, t₇₂=4.9 [P<.001], 2.4 [P=.02], 2.7 [P=.008], 2.7 [P=.009], 2.6 [P=.01], 4.7 [P<.001], 3.7 [P<.001], visit 2: CH 2, t₇₂=2.0 [P=.047], CH 6-7, t₇₂=3.2 [P=.002], 2.6 [P=.01]; VRB group, visit 1: CH 2, t₇₂=2.1 [P=.04], CH 7, t₇₂=4.6 [P<.001]).

Next, for the TMB group, we found activations in the PMC: visit 1, CH 13, t₇₂=2.4 (P=.02), visit 2: CH 11, t₇₂=2.2 (P=.03), CH 13, t₇₂=3.0 (P=.004), CH 14, t₇₂=3.4 (P=.001), CH 26, t₇₂=2.6 (P=.01), CH 33, t₇₂=2.3 (P=.03). We also observed deactivations in the bilateral auditory cortices (visit 1: CH 22, t₇₂=−3.7 [P<.001], CH 23, t₇₂=−2.3 [P=.02]; visit 2: CH 22, t₇₂=−3.2 [P=.002], CH 23, t₇₂=−3.6 [P<.001], CH 42, t₇₂=−4.0 [P<.001]) and right temporal-parietal junction (TPJ; visit 1: CH 39, t₇₂=−2.9 [P=.004]). For the VRB group, we observed activations in the PMC/SMA (visit 1: CH 9, t₇₂=2.8 [P=.007], CH 11, t₇₂=2.8 [P=.007], CH 33, t₇₂=4.3 [P<.001], CH 34, t₇₂=3.6 [P<.001], visit 2: CH 10, t₇₂=2.3 [P=.03], right DLPFC, visit 1: CH 32, t₇₂=7.0 [P<.001], visit 2, CH 32, t₇₂=4.6 [P<.001]) and V1 (visit 2 CH 43, t₇₂=2.2 [P=.03]).

Finally, we performed a contrast analysis between the two groups by combining the two visits, as shown in Figure 7. The results suggested that the VRB group showed greater activation in the right DLPFC (CH 32, t₇₂=5.0 [P<.001, false discovery rate (FDR) corrected]), right TPJ (CH 39; t₇₂=3.4 [P=.001, FDR corrected]), and left STG regions (CH 23, t₇₂=2.7 [P=.048, FDR corrected]). A detailed list of brain activation during the tQST sessions is provided in Multimedia Appendix 2.

It is commonly assumed that MB exerts its analgesic effect through abstract interoception [5]. Interoception refers to consciously refocusing the mind’s attention to the physical sensation of organ function. On the other hand, VR provides participants with an immersive visual-auditory sensory-exteroceptive experience that modulates pain. In this study, we dissected these central analgesic processes by imaging the brains of 2 groups of healthy subjects using fNIRS, exposed to either TMB or a VRB protocol.

Our first finding was that both groups attained a raised pain threshold after one week of breathing practice without significantly different HRV measurements. The results suggested that both TMB and VRB techniques effectively increased pain thresholds in the participants. It is worth noting that the pain-evaluation (tQST) sessions were conducted after the breathing practices, which reproduced their application in the clinical environment, as before a medical/dental procedure. Thus, the analgesic effects were postbreathing effects. The 0.6 °C pain threshold increase in both groups is considered significant, as the human nociceptive sensation thresholds for warmth can be as low as 1.5 °C above the baseline at around 30 °C [38,39]. Nevertheless, the collected PANAS scores suggested that the analgesic mechanisms of the two breathing practices were different. We found that the serenity score increased in the TMB group, although it was accompanied by an increased fatigue score. In contrast, the VRB group showed no significant change in the serenity scores. Serenity is a mental state of being calm, peaceful, and untroubled [40]. This mental state has been shown to increase after meditation [41] and reduce pain [42].

Following the observed temperature threshold changes, we first analyzed the brain activation patterns during the tQST sessions. Our results first confirmed the activation of the contralateral S1 region to noxious thermal stimulation. The S1 has been studied intensively for its critical function in pain and intensity processing [43,44]. As expected, we observed consistent activations in the contralateral (right) S1 region during the tQST sessions in both groups, evoked by the noxious heat stimulation applied to the participants. Interestingly, we noticed a qualitative association between aPFC activation and average pain thresholds across visits and groups, as indicated in Figure 6. Specifically, the TMB group had the lowest pain threshold of 45.4 °C on day 1, accompanied by significant activation in seven out of eight data CHs, whereas the VRB group had the highest pain threshold of 47.1 °C on day 7, with aPFC activation. As one of the cortical executive regions, the aPFC plays an essential role in pain appraisal [45]. Structurally, previous studies reported that chronic pain led to gray matter loss [46-49], whereas meditation increased gray matter volume in the aPFC region [50]. Functionally, studies have found greater aPFC activation with meditation practice [51]. In the context of pain, aPFC activation was correlated with the unpleasantness aspect of pain [52]. However, during the pain process, fMRI imaging revealed reduced functional brain activation in the aPFC region [11,53]. In this study, the trend of less activation was associated with an increased pain threshold, suggesting that less appraisal of pain was induced by the week-long breathing practice. We also found right TPJ, right DLPFC, and V1 activation in the VRB group, whereas TPJ and STG (auditory) deactivation were observed in the TMB group. The TPJ region serves as a hub for integrating multisensory body-related information, including touch, visual, and auditory inputs [54]. In this study, different from the proposed distraction mechanisms [16,18], we observed a postbreathing effect, in which the participant’s TPJ and visual activation were present under the pain condition, even without the VR experience. The right DLPFC possibly provided a modulatory effect on pain to achieve a higher pain threshold. The DLPFC is a critical region that directs attention away from pain [55,56] and inhibits both affective and sensory aspects of pain in the brain [57], as revealed by previous fMRI studies.

To further study the brain mechanisms during the two types of breathing practices, we then analyzed their associated functional connectivity. The TMB group demonstrated a closer working relationship among areas, including the aPFC, PMC, and S1 regions. In contrast, we observed fewer connections to the visual and auditory regions in the VRB group. The immersive 3D, sensory-exteroceptive, virtual experience reinforced the participants’ cortical audio-visual activations, thus depriving the S1 processing of the ascending pain inputs. To study whether functional connectivity is related to pain sensitivity, we examined the correlations between the observed connectivity strength and the pain thresholds measured during the tQST session. Interestingly, we found that the connections among the aPFC, PMC, and S1 regions in the TMB group were associated with pain thresholds across visits. In contrast, in the VRB group, we found that the functional connections among the auditory/visual regions, PMC, inferior parietal lobe, and S1 were associated with temperature thresholds. A previous resting-state fMRI study found lower pain sensitivity in meditators with decoupled executive and pain-relevant brain regions [11]. Similarly, in our study, we found lower pain sensitivity in the TMB group participants with less connectivity between the aPFC and PMC-S1 regions. For the VRB group, lower pain sensitivity was observed in participants with decoupled visual-auditory-DLPFC and PMC regions.

On the basis of these findings, we propose two possible mechanisms for the TMB and VRB practice—in the TMB group, the aPFC modulated attention [58-60] and contextual evaluation of internal sensory events [61-63]. The PMC, as part of the mirror neuron system [64], increased its functional connection with the S1 to facilitate the sensory-interoceptive processing of breathing. This process inhibited the S1 in sensory-discriminative pain processing during later tQST sessions. In contrast, VR induced an immersive 3D exteroception with augmented visual-auditory cortical activations to diminish functional connection with the S1, consequently weakening the pain processing function of the S1.

There were some caveats in this study. First, we asked the participants to mimic the VR experience for at-home practices between day 1 and day 7 lab visits. However, the at-home practice was not as immersive as the in-lab VR practice (owing to the noncompletely portable apparatus of the technologies), which might have dampened its effect. We will use a mobile VR device that works with smartphones in our future studies to address this issue. Next, instead of a nonintervention control group, we used the TMB group as the active control group. Although we found different brain activation and connectivity patterns, in contrast to the VRB group, we still need to compare both groups with a nonintervention control group and evaluate the effects of both the TMB and VRB breathing practices in the future.

In conclusion, as shown in Figure 8, our study suggested two distinct analgesic mechanisms of VRB and TMB practices, following the concept of dualism, Yin and Yang, in ancient Asian philosophy [65]. On the Yin-side, the aPFC, activated by the TMB practice, modulated the PMC to maintain an uninterrupted sensory-interoception via the S1, which prevailed over its sensory-discriminative processing of the ascending pain. On the Yang-side, the VRB practice brought an immersive 3D sensory-exteroceptive VR experience via augmentation of cortical visual-auditory activations that overrode the pain processing function to raise the pain threshold.