Decoding spatial location of perceived pain to acupuncture needle using multivoxel pattern analysis

Abstract

The present study applied multivoxel pattern analysis to decode spatial discrimination in pain perception to acupuncture needle from brain functional magnetic resonance image. Fourteen participants were stimulated by acupuncture needles at two adjacent body parts on their left forearm (PC6 vs. HT7). We trained support vector machines on the spatial information from the whole-brain functional magnetic resonance imaging data and projected the support vector machine weight to the brain image space to represent the effect of each voxel on the classifier output. Using region-of-interest masks in individual brains, we trained and tested a linear support vector machine classifier on the accuracy of spatial discrimination in trial-wise functional magnetic resonance imaging data. A classical univariate general linear model analysis testing for differences between the two different locations did not reveal any significant differences. Multivoxel pattern analysis revealed that the brain regions for the prediction of sensory discrimination in pain perceptions to two different points were in the primary somatosensory cortex, primary motor cortex, and supramarginal gyrus, anterior and posterior insula, anterior and posterior cingulate cortex, ventromedial prefrontal cortex, and inferior parietal lobule. Our findings suggest that spatial localizations of pain perceptions to acupuncture needle can be predicted by the neural response patterns in the somatosensory areas and the frontoparietal areas.

Introduction

The central mechanisms associated with somatic localization of tactile stimuli are well established in the field of neuroscience.^1–6 The primary somatosensory cortex (SI) is a key sensory receptive area for somatic stimuli, ⁷ and perceptual discrimination of somatosensory stimulation is encoded in its neurons. ⁸ The parietal cortex is also involved in the neural processing of information concerning the location and intensity of somatosensory stimuli. ⁹ A functional magnetic resonance image (fMRI) study on a two-point discrimination task revealed that cognitive discrimination of spatially distinct stimuli is associated with the supramarginal gyrus (SMG) of the inferior parietal lobule (IPL). ¹⁰ Recently, using advanced machine learning techniques, multivoxel pattern analysis (MVPA) of fMRI data successfully showed that brain activation in the SI contains more distinctive spatial patterns that decode the location of tactile stimuli placed closely together on the skin surface, whereas secondary somatosensory cortex (SII) and SI have equal accuracy in decoding widely spaced tactile stimuli. ² The hierarchical view of a somatosensory system ranging from unimodal somatosensory functions to higher order cognitive brain areas can be applied to understanding the process of somatosensory localization on the body scheme.

Somatosensation includes the processing of tactile, proprioceptive, and nociceptive information. ¹¹ Both tactile and painful stimuli produce brain activation in similar regions of the SI and SII. Painful stimuli are also associated with the anterior insula and frontal cortices, or regions highly linked to the limbic systems and emotional processing.^12–14 The distinct modules are differentially engaged in discrimination of sensory features of nociceptive information. ¹⁴ Much like other sensory modalities, intensity-related information of pain perception is preferentially processed by ventrally directed processing stream, while spatial information of pain perception is preferentially processed by the dorsally directed processing stream.^14,15 The somatotopic representations of nociceptive information are useful in explaining the behavioral effects of spatially directed placebo analgesia. ¹⁶ The site of nociceptive stimulation, for example, arm or leg, could be decoded from brain patterns during the anticipation and perception of painful stimulation. ¹⁷

Acupuncture action is known to exert by a site-specific action.^18,19 The issue of point specificity has been one of the critical issues in the field of acupuncture research. ²⁰ Previous studies applied general linear model (GLM) approach to investigate neural response to acupuncture stimulation using a set of predefined regressors. On the other hand, MVPA method can be expected to overcome a lack of sensitivity of mass-univariate approaches, computing voxelwise statistics that have been applied in most studies of acupuncture.^17,21 Since acupuncture needle stimulations can be applied to two adjacent sites, it would be more important to distinguish the brain patterns to needle stimulations at different acupoints. Despite the importance of understanding discrete spatial information of the source of pain, the central mechanism that allows spatial discrimination of painful sensations is not fully understood.

Here, two adjacent body parts on the left forearm (median vs. ulnar nerve) were stimulated by acupuncture needle during fMRI scanning. The MVPA methods could provide considerable increases in the amount of information than the traditional univariate methods. ²² Prior to the multivariate analyses, we performed univariate analyses for comparison purposes. We applied MVPA to decode spatial discrimination of painful stimulations on two discrete locations on the forearm from human brain signals.

Methods

Participants

Fourteen right-handed male participants (mean age 22.1 ± 1.1 years) took part in the study. They had no known history of neurological, psychiatric or visual disorders. Participants were prohibited from drinking alcohol or caffeine and from taking any drugs or medications on the day of the experiment. After informing them the nature of the experiment, they provided full written consent. The study was conducted in accordance with the Declaration of Helsinki and was approved by the Institutional Review Board at Korea University.

Experimental design

Using a multivariate method, we revisited a data set from a previous study. ²³ The previous study dealt with commonalities and differences in brain responses to enhanced bodily attention around acupuncture points with and without stimulation. ²³ In the present study, however, we used MVPA and revealed brain activity patterns encoding spatial discrimination in the forearm in pain perceptions. MVPA has increased sensitivity for the detection of cognitive states, compared with the univariate method. ²⁴ On the other hand, the target of the previous study was not the location of painful stimulation. Therefore, points A and B were not distinguished from each other, but rather grouped as stimulation in the analysis. In addition, the previous study had two sessions: the first session provided actual stimulation, while the second session induced attention without actual stimulation. However, this study employed only the data of the first session. Here, we provide a summary of the experimental design, focusing on aspects relevant to the question being addressed in the present study (Figure 1).

Figure 1.

Experimental design. Twenty trials were conducted during the experiment. Location A (median nerve) and location B (ulnar nerve) were located in the forearm and were separately innervated in each position. In a trial, they were randomly selected and stimulated by rotating an acupuncture needle. Locations A and B were stimulated 10 times. In a trial, stimuli were given for 6 s after a resting period of 16 s, and then 8 s were given to report the response to the stimuli, including location and intensity. Based on the regression model for each trial, the coefficient maps were used as samples for the SVM classifier to be trained and tested. SVM: support vector machine.

During the fMRI scanning, participants were asked to focus on the stimulated area on their left forearm. The whole session of the experiment consisted of 20 trials. Each discrimination trial started with a 16-s rest period, during which participants were told to fixate on a red cross, followed by a 6-s acupuncture stimulation period at either location A or B (PC6: median nerve innervation; HT7: ulnar nerve innervation, respectively). During this period, a blue fixation cross was shown as the visual stimulus. The order of locations within each session was counter-balanced and randomized among the participants. Following the stimulus, participants were required to identify the stimulated locations by pressing one of two buttons representing points A and B on a four-button MRI-compatible button-box (current design) and then to rate the intensity of the sensation by pressing one of the same button-box, which was held in the right hand during the session. Each rating screen appeared for 4 s, and participants were instructed to make their decision before the end of this period. The overall time duration of each trial was 30 s. The performance data of the discrimination task by participant were included in our previous study. ²³ Among 14 participants, 11 participants were 100% accurate, and the other 3 participants were 90%–95% accurate. There were only four wrong responses. Since this study aimed to decode the physically stimulated location, we did not exclude wrong responses from the data samples.

For the fMRI acquisition series, participants’ discrimination responses were recorded using Matlab Psychtoolbox. Data were processed using customized programs within the R software package. The detailed design is fully described in our previous paper. ²³

Data acquisition and preprocessing

The fMRI scans were acquired using a 3T Siemens Tim Trio magnetic resonance device, with a head coil attached. To minimize movement artifacts, the head of each subject was stabilized using a head holder, and all images were acquired by a well-trained professional operator. In each scan session, 300 volumes of the entire brain were collected in 37 axial slices (repetition time (TR) = 2000 ms, echo time (TE) = 30 ms, flip angle = 90°, field of view = 240 × 240 mm², voxel size = 3.8 × 3.8 × 4.0 mm³). As an anatomical reference, a three-dimensional T1-weighted magnetization-prepared rapid gradient echo image data set was acquired using the following parameters: TR = 2000 ms, TE = 2.37 ms, flip angle = 9°, field of view = 240 × 240 mm², voxel size = 0.9 × 0.9 × 1.0 mm³, and 192 slices. The Echo Planar Imaging (EPI) images were corrected for slice timing and realigned to the first volume using sinc interpolation. EPI images were coregistered to the structural T1 images. The images were transformed to a common space (Talairach space). ²⁵ The spatial smoothing process was applied separately between univariate analysis and multivariate analysis, and no spatial smoothing was applied in multivariate analysis. Images showing motion or rotation greater than 3 mm or 3° were excluded. Excluded images were less than two percentage of the whole data set.

Univariate analysis

Prior to the multivariate analyses, we performed univariate analyses for comparison purposes. For the univariate analyses, surface-based analysis was applied. FreeSurfer software was used to separate each anatomical volume into gray and white matter structures and to perform cortical surface reconstruction. The cortical surfaces in the EPI images were smoothed using a Gaussian filter with a full width at half maximum of 6 mm. For each location-specific pain perception (locations A and B), a boxcar function was used to represent the duration of each event and convolved with a gamma function. Contrast images corresponding to stimulation or responses to “location A” and “location B” were generated by fitting “location A” and “location B” regressors to scan time courses using the analysis of functional neuroimages (AFNI) program 3dDeconvolve. ²⁶ A “summary statistics” procedure involving one-sample t-tests performed across individual contrast images was used to model group effects. Cluster threshold criteria were determined using Monte Carlo simulations, which resulted in a family-wise error (FWE)-corrected significance threshold of p < 0.05 ²⁷ using AFNI AlphaSim program (http://afni.nih.gov/afni/docpdf/AlphaSim) with voxel-wise statistical threshold (t-score > 2; p < 0.032).

Whole-brain classifier weight analysis

To identify the brain regions responsible for discriminating two locations, we examined support vector machine (SVM) weight values on each voxel from the whole-brain SVM classifier. The trial-specific estimates were obtained through a beta-series regression, which constructs an estimation model with separate regressors. ²⁸ The trial-wise regressors were modeled using gamma hemodynamic response function. Separate parameter estimate images (beta images) for each trial were obtained through these processes, and parameter estimates were obtained for each participant. The resulting parameter estimates were processed further for the multivariate analysis using Nibabel (http://nipy.org/nibabel) and Scikit-learn software. ²⁹ First, we normalized trial-specific parameter estimates to achieve centering relative to the mean and unit variance. Second, voxels were selected from a brain mask generated for the corresponding participant using AFNI. To test the performance of the classifier, we used leave-one-trial-out cross validation. Specifically, we selected one trial from the 20 trials in the session to omit from training. One more trial from the class opposite to that of the omitted trial was selected randomly and also excluded from the training.

We explored brain regions containing information associated with somatic discrimination in the pain perception. The SVM algorithm was applied to map whole trial-specific data without prior selection of features. By treating the whole brain as a point in a high-dimensional space, the SVM linearly could classify trial-specific beta images into two classes (perceptions of stimulation at location A and location B) by finding an optimally separating hyperplane, which is determined by applying a linear function separating the training data with maximal margin. The optimal hyperplane is trained with a weight vector that indicates the direction of perceptions from which the two stimulated locations differ. The weight vector represents a pattern of the most discriminating voxels. ³⁰

Apart from the cross validation for the accuracy evaluation, SVM weights were extracted from the classifiers which were trained using all data sets without distinction between training set and validation set. The whole-brain data were used to train the SVM classifier and extract the weights for each subject. Each element of the SVM weight vector corresponds to voxels of the whole brain. The SVM weight vector encodes the contributions of all voxels to the classifier. The absolute size of the SVM weight relative to other voxels gives an indication of how important the feature was for classification. Before evaluating statistical significance of the weight of each voxel, we took the absolute value of the weights.

We aimed to evaluate significantly influencing brain areas on the classification of stimulated location across participants. We conducted two permutation tests: one for accuracy (trained using leave on trial out) and one for weight (using all data). Group mean of each voxel weight was evaluated using one sample t-test against null distribution of weights sampled during the permutation test for the whole-brain classification accuracy. In our study, null distribution of SVM weight in each voxel was estimated as a Student’s t-distribution calculated from the sample mean and sample variance of the group mean weights obtained from the 1000 permutations. The Student’s t-distribution was calculated with N − 1 degrees of freedom where N is the number of participants. ³¹ To obtain null distribution of group mean accuracy, we first obtained all required null accuracies (14,000 null accuracies = 14 participants × 1000 iteration) and calculated 1000 group mean accuracies by averaging across participants. Then, we probed where the “real” weight t-value falls on this null distribution. While other studies employ the average across participants per iteration, we considered that the order of averaging process was of less importance since the aim of this procedure was to generate a distribution under the null hypothesis. We conducted the same procedures for permutation test for the region of interest (ROI) analysis and the weight analysis.

Cluster threshold criteria were determined using Monte Carlo simulations, which resulted in a FWE-corrected significance threshold of p < 0.05 ²⁷ using AFNI AlphaSim program (http://afni.nih.gov/afni/docpdf/AlphaSim) with voxel-wise statistical threshold (t-score > 2; p < 0.017). Estimation of spatial smoothness of the data and followed iterative t-statistic map simulations were processed on 3D whole-brain data not on the surface. The resulted statistical map in Talairach space was projected to a surface (3dVol2Surf) and visualized in an inflated surface using AFNI SUMA.

ROI classification analysis

To assess and compare the degree of contribution of the brain regions to the process of encoding the localization of painful stimuli, further analysis was performed using ROIs. The ROIs used for MVPA were brain regions parcellated according to Desikan–Killiany–Tourville protocol.^32,33 These regions were defined on the basis of anatomical landmarks, independently from the results of the above whole-brain analysis. Desikan–Killiany–Tourville protocol distinguishes cortex into 32 areas for each hemisphere. ³³ Automatically assigned neuroanatomical labels for each location on the cortical surface were used to define ROIs. The reconstructed surface of each participant was resampled from FreeSurfer to AFNI SUMA’s standard mesh topology using MapIcosehedron (200,000 triangles; 100,002 nodes), and voxels corresponding to FreeSurfer’s cortical parcellation (Desikan–Killiany–Tourville Atlas) were extracted using ROI masking.^32,33 In contrast to the first multivariate analysis using whole-brain images, we applied the SVM algorithm to each ROI in the second analysis. We evaluated performance accuracy using the same leave-one-trial-out cross validation paradigm used in the whole-brain analysis. The process was repeated for every participant, and mean classification accuracies were calculated among participants for each ROI.

We used a permutation test to examine whether or not classification accuracy exceeds the accuracy attributed to chance (50%). The statistical significance of the classification accuracy from each ROI was evaluated through nonparametric permutation test. The acupuncture stimulation labels were randomly shuffled for each ROI. Predicting stimulated labels using SVM classifier trained with the permutated data sets was repeated 1000 times. The statistical significance of group mean accuracy obtained from each ROI was evaluated by comparing with null distribution of group mean accuracies accumulated from 1000 permutations. The null distribution of group mean accuracies for each ROI was estimated in the same way with whole-brain classifier. Since multiple ROIs were analyzed, multiple comparisons correction with false discovery rate <0.05 was performed using the Benjamini–Hochberg method.

Results

Behavioral results

In pain perception, there were no significant differences between location A (2.50 ± 0.89) and location B (2.65 ± 0.94) (t = 0.701, p = 0.496). The other behavioral results were of no interest in the present analysis (see a previous study for a full description ²³ ). The correct rates were up to 98.6%.

Univariate analysis

A classical univariate GLM analysis testing for differences between the two different locations (location A vs. location B) did not reveal any significant differences at a threshold of p < 0.001 (uncorrected for multiple comparisons). Univariate analysis revealed that stimulation of both the location A and location B produced brain activation in the bilateral insula, operculum, inferior frontal gyrus, supplementary motor area (SMA), SI and SMG, and deactivation in the default mode network (DMN) consisting of the ventromedial prefrontal cortex (vmPFC), PCC, IPL, medial temporal gyrus (MTG), and parahippocampus at an FWE-corrected significance threshold of p < 0.05 (Figure 2).

Figure 2.

Univariate analysis. In the univariate analysis, stimulation of both locations (location A and location B) produced brain activations in the bilateral insula, operculum, inferior frontal gyrus, supplementary motor area, primary somatosensory cortex and supramarginal gyrus and deactivations in the default mode network, consisting of the ventromedial prefrontal cortex, posterior cingulate cortex, inferior parietal lobe, medial temporal gyrus, and parahippocampus. FWE: family-wise error.

Whole-brain classifier weight analysis

To map the discriminative information for spatial locations stimulated by noxious stimuli (a rotating acupunctured needle), we trained and tested a linear SVM classifier on trial-by-trial correlates of whole-brain fMRI data when location A or B was stimulated. We found that monitoring whole-brain activity during the pain perception (58.6%) enabled statistically significant predictions (p < 0.001; nonparametric permutation test using N = 1000 permutations). One sample t-test of SVM weight map against null distribution sampled from permutation test revealed that SI, primary motor cortex (MI), paracentral cortex, anterior and posterior insula, SMG, anterior cingulate cortex (ACC), vmPFC, PPC, and IPL were significantly important regions that allowed statistically significant predictions of the two stimulated sites at a FWE-corrected significance threshold of p < 0.05 (Figure 3).

Figure 3.

Whole-brain classifier weight analysis. In the multivariate analysis, the SI, MI, paracentral cortex, anterior and posterior insula, SMG, ACC, vmPFC, PPC, and IPL allowed for statistically significant discrimination between the two stimulated sites. The percentages on the right indicate the resulting classification accuracies for each session (p < 0.001). FWE: family-wise error; MI: primary motor cortex; SMG: supramarginal gyrus; SI: primary somatosensory cortex; ACC: anterior cingulate cortex; vmPFC: ventromedial prefrontal cortex; PPC: posterior parietal cortex; IPL: inferior parietal lobe.

ROI classification analysis

All cortical ROIs were selected based on Desikan–Killiany–Tourville protocol.^32,33 Our analysis generated a ranked order of brain regions showing significantly higher accuracy than chance-level (50%) on spatial discrimination. Among them, the MI (contralateral to stimulation site; 65%), SMA (contralateral to stimulation site; 64%), SMG (contralateral to stimulation site; 62%), SI (contralateral to stimulation site; 62%), and dorsolateral prefrontal cortex (dlPFC; contralateral to stimulation site; 62%) showed accurate trial-by-trial discrimination in terms of pain perception (Figure 4).

Figure 4.

ROI classification analysis. Brain regions which have significantly higher classification accuracies are shown with null simulated distribution. The MI (contralateral to stimulation site; 0.646), SMA (contralateral to stimulation site; 0.643), SMG (contralateral to stimulation site; 0.621), SI (contralateral to stimulation site; 0.621), and dlPFC (contralateral to stimulation site; 0.621) showed accurate trial-by-trial discrimination in terms of pain perception. dlPFC: dorsolateral prefrontal cortex; MI: primary motor cortex; SMA: supplementary motor area; SMG: supramarginal gyrus; SI: primary somatosensory cortex.

Discussion

The present study statistically assessed each set of multivoxel patterns in terms of the perception of spatial information of pain in the two adjacent body parts. Our MVPA-based approach in the current study revealed distinguishable brain activation patterns during spatial discrimination in pain perception with significantly higher accuracy than chance level. According to the classification performances for each brain region, the somatosensory processing regions, such as the SI, and frontoparietal brain areas, including the SMG and dlPFC, were main regions featuring distinctive brain activation patterns in the discrimination of two different stimulated locations.

Our study applied MVPA methods and revealed that spatial information of pain perception could be decoded from neural activation patterns in brain areas including SI, dlPFC, and SMG. These neural patterns in the current study were consistent with the pain perception suggested by the following previous studies. Localization of somatosensory stimulation extending far beyond SI involves higher level somatosensory processing areas including the SII, SMG, IPL, and PPC in the hierarchical processing scheme.^1,2,5,6 Similar to the processing of visual and auditory, and other innocuous somatosensory information, frontoparietal interactions were critically involved in the discrimination of spatial features of pain. ¹⁵ The two different sensory features of noxious stimuli (intensity and location) are preferentially processed by distinct neural system. Ventrally located modules, such as insula cortex and cingulate cortex, are much associated with intensity-related information of pain perception, whereas a dorsally located module, just as PPC and right dlPFC, is much associated with spatial information of pain perception. ¹⁴ In this study, frontoparietal areas, including SMG and dlPFC, were also preferentially prominent during spatial discrimination of needles at two adjacent acupoints. Frontoparietal areas beyond the somatosensory areas contribute to spatial discrimination of noxious needle stimuli. Ritter et al. reported that the spatial information associated with painful stimuli was decodable from brain patterns of the prefrontal cortex, rostral ACC, and parietal operculum. ¹⁷ Although the two discriminated locations in the current study were very close to each other on the forearm (median vs. ulnar nerve stimulation) compared with other previous studies (arm or leg), nociceptive location could be predicted successfully from a broad network of sensorimotor processing regions.^15,17

The somatosensory system is composed of hierarchical structures such as thalamus, SI, and SII, and it is known that discriminative somatosensory information is distributed in various brain regions. ⁶ The somatosensory processing regions, which are responsible for the localization of nonpainful tactile stimulation in previous studies, were common to brain regions involved in the localization of acupuncture stimulation in this study. Motor areas including MI and SMA were detected for the spatial information regarding the site of needle stimulation. Although no actual movements were needed during needle stimulation, and discrimination-related brain patterns were analyzed before subjects’ response, someone might raise the concerns about which motor areas’ activations might be associated with motor-related processing, such as motor planning for pushing the switch to indicate a choice. However, most of brain regions that we found from ROI analysis were more dominant in contralateral to the stimulation sites (i.e., ipsilateral to the response site). Although ipsilateral signals also reflect an efference copy of a contralateral motor command and facilitate coordination between both limbs in bimanual behavior, brain patterns within frontoparietal areas are mainly implicated in contralateral movement planning.^22,34 Thus, it is more likely that different brain activity patterns in our study may be related to the spatial discrimination of pain perception.

The whole-brain accuracy, including all the information of these regions, was not greater than the accuracy of the individual ROIs in the current study. One of the reasons seems to be due to the limitation known as “curse of dimensionality.” The voxel number of whole-brain data is estimated about 35,000. When the dimensionality of feature space increases, the volume of the space increases so fast that the useful data might become sparse. The other reason might be derived from redundant information of the brain regions included in the whole-brain analysis, which do not have relevant activities to decode the painful stimulation site. To obtain higher accuracy of whole-brain data in this study, therefore, additional feature extraction method should be considered in the future studies.

The traditional univariate GLM methods are known for insufficient sensitivity for more complex processes such as discrimination of stimulus locations, which often involves widely distributed neuronal activities. In contrast, multivariate analysis could detect different brain patterns to stimulus parameters and accumulate the weak information available at each brain area in an efficient way. ³⁵ Similarly, we could not identify distinct brain regions to two adjacent body parts (median vs. ulnar nerve; PC6 vs. HT7 acupoints) by using the conventional GLM methods. Using the more sensitive approach of the MVPA method, we were successful in decoding brain regions that showed predictable activity patterns in spatial information of pain perception to acupuncture needle stimulation of two different sites. A previous study rigorously employed the conventional univariate GLM methods to demonstrate that acupuncture at vision-related acupoints (such as BL60 and GB37) specifically produced brain activation in the occipital cortex, of which the results however were unsuccessful. ³⁶ It is believed that the brain activation and deactivation patterns in response to acupuncture needle stimulation reflect the sensory, cognitive, and affective dimensions of pain. ²⁰ Conversely, Li et al. applied MVPA and was able to classify brain activity patterns that differed between a vision-related acupoint and a nearby nonacupoint. ²¹ The debate concerning the existence of acupoint specificity continues to be solved in future research.

Brain activations to acupuncture stimulations reflect not only the pain perception to the stimuli but also various cognitive and emotional responses to the stimuli, expectations to the treatment, and physiological actions of the treatment. The previous studies demonstrated that acupuncture stimulations bilaterally produced activations in sensorimotor cortical network such as insula and ACC and deactivations in DMN, such as medial prefrontal cortex and parahippocampus.^20,23 In the current study, a classical univariate GLM analysis showed bilateral activations in the sensorimotor cortex and deactivations in the DMN in both of the acupuncture stimulated locations (Figure 2). On the other hand, ROI classification analysis revealed that spatial discriminations are mainly involved in the contralateral frontoparietal brain regions (Figure 4). Since decoding accuracy represents the degree to which a condition is represented in the pattern of activity distributed across a region’s multiple voxels, the MVPA can provide sensitivity to information that cannot be detected in mean activation levels alone. ³⁷ A recent study demonstrated that the multivariate measure successfully detected the lateralization of orthographic processing in the visual word form area. ³⁸ It is assumed that our analysis on the spatial discriminations of pain using MVPA reflects lateralization of brain function, that is, predominant in the contralateral side than GLM analysis. Speculatively, we propose that multivariate techniques can be more useful to understand how the lateralized and bilateral brain achieves pain perception. However, further work is required which relates other behavior measures to multivariate laterality.

Several limitations should be noted in this study. Although there were no significant differences of subjective pain ratings between the two locations in the forearm, we could not rule out the possibilities of the involvement of different intensities of pain perception to decode spatial discrimination of painful stimulations from the two different points. The classifier implemented in this study was purposed to predict stimulated location by acupuncture. No control condition was required to reveal whether a brain region contains a neural response discriminating the stimulated location. However, comparing with tactile stimuli would benefit to find more specific brain regions for acupuncture stimulation. We did not include the tactile stimuli as a control. Thus, it will be necessary to further explore the prediction performance of spatial patterns of brain activity during tactile stimulation in the future study. The somatosensations from acupuncture needling include a constellation of sensations experienced by patients such as heaviness, soreness, numbness, and distension and may have affective connotations such as feeling refreshed or relieved.^39,40 The main manifestations of deqi sensations are in general separated from the acute pain at the site of the needling, especially in the case of sharp pain. ⁴¹ Since we did not identify the characteristics of deqi sensations from the participants’ responses in this study, we were not able to clearly differentiate the brain patterns of pain perceptions from other kinds of somatosensations to acupuncture needle. Further study is necessary to characterize the complex sensations of acupuncture needling. Furthermore, MVPA analyses in the current study were applied using only activation patterns estimated from single trials within the same run. Inflated false positives can be minimized when the testing and training sets in cross validation do not contain patterns estimated from the same run. ⁴² Last but not least, we used a less conservative t value for voxel-wise statistical threshold and it might be susceptible to false positives. ⁴³

In summary, our results show that the site of nociceptive information from the two nerve innervations in the forearm can be successfully decoded from spatial patterns of brain activity during needle stimulation. Our findings suggest that spatial information of pain perceptions to acupuncture needle is represented in the somatosensory processing regions as well as frontoparietal brain areas, such as SMG and dlPFC. We strongly believe that these findings could offer new insights into understanding brain processing of spatial information contained in human pain perception.