NeuroImage 90 (2014) 153–162

Contents lists available at ScienceDirect

NeuroImage journal homepage: www.elsevier.com/locate/ynimg

Ranking brain areas encoding the perceived level of pain from fMRI data Stefania Favilla a, Alexa Huber a, Giuseppe Pagnoni a, Fausta Lui a, Patrizia Facchin a, Marina Cocchi b, Patrizia Baraldi a, Carlo Adolfo Porro a,⁎ a b

Department of Biomedical, Metabolic and Neural Sciences, University of Modena and Reggio Emilia, Via G. Campi 287, Modena, Italy Department of Chemical and Geological Sciences, University of Modena and Reggio Emilia, Via G. Campi 183, Modena, Italy

a r t i c l e

i n f o

Article history: Accepted 1 January 2014 Available online 10 January 2014 Keywords: Cerebral cortex Multivariate analyses Partial least squares regression Pain perception Thalamus Variable importance of projection

a b s t r a c t Pain perception is thought to emerge from the integrated activity of a distributed brain system, but the relative contribution of the different network nodes is still incompletely understood. In the present functional magnetic resonance imaging (fMRI) study, we aimed to identify the more relevant brain regions to explain the time profile of the perceived pain intensity in healthy volunteers, during noxious chemical stimulation (ascorbic acid injection) of the left hand. To this end, we performed multi-way partial least squares regression of fMRI data from twenty-two a-priori defined brain regions of interest (ROI) in each hemisphere, to build a model that could efficiently reproduce the psychophysical pain profiles in the same individuals; moreover, we applied a novel three-way extension of the variable importance in projection (VIP) method to summarize each ROI contribution to the model. Brain regions showing the highest VIP scores included the bilateral mid-cingulate, anterior and posterior insular, and parietal operculum cortices, the contralateral paracentral lobule, bilateral putamen and ipsilateral medial thalamus. Most of these regions, with the exception of medial thalamus, were also identified by a statistical analysis on mean ROI beta values estimated using the time course of the psychophysical rating as a regressor at the voxel level. Our results provide the first rank-ordering of brain regions involved in coding the perceived level of pain. These findings in a model of acute prolonged pain confirm and extend previous data, suggesting that a bilateral array of cortical areas and subcortical structures is involved in pain perception. © 2014 Elsevier Inc. All rights reserved.

Introduction Understanding the cerebral mechanisms of pain perception is a challenging issue, with highly relevant theoretical and clinical implications. Acute pain is associated with widespread changes in the activity of brain structures in experimental animals and in humans, as suggested by different kinds of mapping techniques (Apkarian et al., 2005; Duerden and Albanese, 2013; Porro, 2003; Tracey, 2008). These functional changes are conceivably related not only to the perceived aspects of pain, but also to motor and vegetative reactions to acute noxious stimulation (Piche et al., 2010; Vogt, 2005), and to arousal and stimulus saliency (Legrain et al., 2011). Current evidence suggests that pain emerges from the activity of widespread cerebral networks (Ohara et al., 2008; Ploner et al., 2010), among which subsets of regions can be involved in different aspects of pain experience (e.g. sensory–discriminative vs. cognitive– evaluative) (Kulkarni et al., 2005; Oshiro et al., 2009; Wiech et al., 2008).

⁎ Corresponding author at: Department of Biomedical, Metabolic and Neural Sciences, University of Modena and Reggio Emilia, Via G. Campi 287, I-41125 Modena, Italy. Fax: + 39 059 2055363. E-mail address: [email protected] (C.A. Porro). 1053-8119/$ – see front matter © 2014 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.neuroimage.2014.01.001

However, identifying the most relevant brain regions for pain perception poses various methodological problems. On the one hand, disentangling perceptual brain correlates from changes in brain activity related to vegetative, arousal or motor reactions is more difficult when using brief noxious stimuli, characterized by high saliency, which have been adopted in many experimental pain studies. On the other hand, brain mapping techniques such as functional magnetic resonance imaging (fMRI) yield massive multivariate datasets, from which relevant information is only partially extracted using voxel-based conventional analyses. Different statistical methods have been recently applied, and proven to be useful for decoding and predicting pain from brain activity (Brodersen et al., 2012; Marquand et al., 2010; Prato et al., 2011; Schulz et al., 2012; Wager et al., 2013). Multi-way partial least squares regression (N-PLS) (Bro, 1996) is an extension of partial least square regression (PLS), extensively used in chemistry (Pardo et al., 2004; Pravdova et al., 2002), and firstly applied to functional neuroimage analysis by McIntosh et al. (1996) with the goal of extracting commonalities between brain activity and behavior or experimental design. These regression tools belong to latent variables based techniques and could be used to identify robust patterns of covarying neural activity, and to relate these patterns to task behavior. The usefulness of these approaches in the interpretation of fMRI data

154

S. Favilla et al. / NeuroImage 90 (2014) 153–162

has been recently reviewed (Krishnan et al., 2011; McIntosh and Misic, 2013). In this study, we have applied N-PLS to the analysis of a BOLD-fMRI experiment based on a pain paradigm of acute prolonged pain in healthy volunteers, induced by ascorbic acid injection (Rossi and Decchi, 1997). Our goal was to build a model that could extract the main variation in our data (the time series of fMRI signal changes in brain regions in different volunteers) in the simplest possible way, in order to efficiently reproduce the psychophysical pain profile in the same individual. Furthermore, in order to identify the brain regions mainly involved in determining the perceived pain profile, we have applied a three-way extension (Favilla et al., 2013) of a method, known in the literature as variable importance in projection (VIP) (Wold et al., 1993), useful to summarize the contribution of each variable to the model. The results of this novel approach are compared with those of conventional whole-brain, voxel-based multiple regression analysis, and with the results of a region of interest (ROI)-based analysis on mean estimated values for the psychophysical pain rating regressor. Materials and methods Experimental procedures Seventeen healthy right-handed volunteers (5 males and 12 females, age range 18–40, mean 22.7 years) participated in the study, which was approved by the local Committee on Ethics, after giving their written informed consent. Each volunteer was subjected to a subcutaneous (s.c., n = 8) or intramuscular (i.m., n = 9) injection of an ascorbic acid solution (0.5 ml, 20%), through a plastic cannula inserted (using a 25 gauge needle) into the s.c. or i.m. tissue of the left thenar eminence 10 min before the beginning of the fMRI experiment; therefore, no simultaneous cutaneous stimulus, either tactile or noxious, was given at the onset of the chemical stimulation. The volunteers were instructed to code the sensory intensity of pain throughout the experiment, by moving a computer-controlled horizontal visual analog scale (VAS) using their right (unstimulated) hand. The VAS scale was anchored at “no pain” and “the maximum imaginable intensity of pain”. Raw data from 14 volunteers of this same data set were used in another study by our group, where a completely different approach was used in order to predict the perceived pain intensity in other individuals on the basis of an iterative learning algorithm on fMRI data (Prato et al., 2011). For each volunteer, mean VAS values were computed over consecutive time intervals, each lasting 4 s (same as the TR, see subsection “MRI data acquisition”), and converted into a 0–100 pain intensity scale. In order to construct a suitable predictor for the neural response associated with perceived pain intensity, each subject's psychophysical pain rating curve was also scaled to a 0–1 range and convolved with a gamma function modeling the hemodynamic response. MRI data acquisition The fMRI session lasted 20 min, with the exception of one subject who received the i.m. ascorbic acid injection, for whom the fMRI acquisition lasted only 18 min due to technical reasons. In each subject, a brief “warning” signal (cleaning the skin of the hand with an antiseptic solution for 10 s) was delivered at minute 3 and was followed, 1 min later, by the ascorbic acid injection through the previously inserted cannula. Functional images were acquired by a GE 1.5 T scanner, using an echoplanar (EPI) BOLD sensitive sequence (TR = 4 s; 3.75 × 3.75 × 4 mm voxel size). During the experiment, 300 brain volumes (time points; 270 in one subject) were acquired from 24 contiguous axial planes covering the diencephalon and telencephalon. For each subject, a high resolution T1-weighted image of the whole brain (TR = 35 ms; TE = 9 ms; 120 axial planes: 0.94 × 0.94 × 1.3 mm voxel size), was also acquired for anatomical localization.

ROI selection Twenty-two regions of interest (ROIs) were manually identified bilaterally by the same experimenter (P.F.), and independently checked by another author (C.A.P.), in both hemispheres on the high resolution T1-weighted images of each subject, according to specific landmarks (e.g. Paxinos and Mai, 2004; Vogt, 2005; Yousry et al., 1997). ROIs included 4 subcortical regions (caudate nucleus, putamen, medial and lateral thalamus), and 18 frontal, parietal, cingulate and insular cortical regions: the putative foot, hand, and face representation areas, both of the primary somatosensory (postcentral gyrus) and of the motor (precentral gyrus) cortices, the parietal operculum, the inferior and superior parietal lobules, the precuneus, the medial prefrontal cortex, the anterior, mid- and posterior cingulate cortex, and the anterior and posterior insula. These regions were selected on an a priori hypothesis on their potential involvement in sensori-motor mechanisms and nociceptive processing (Apkarian et al., 2005; Oshiro et al., 2009; Tracey, 2008). Care was taken to include only gray matter in the ROIs. Region inclusion probability maps, showing the across-subject overlap in the MNI space of the 22 anatomical regions of interest (ROIs) of the left hemisphere, are displayed in Supplementary Fig. 1. For each volunteer, the ROI-mask image was then co-registered with the functional volumes, with the following procedure: the T1 image was co-registered to the first volume of the motion-corrected EPI images and the estimated transformation matrix was then applied to the ROI-mask image (using nearest-neighbor interpolation and keeping the original resolution of the ROI mask). The EPI images were resampled to 2 × 2 × 2 mm in order to match the ROI mask. For ROI analyses, we chose to work with the finer ROI resolution instead of the original EPI resolution, in order to avoid the disruption of some of the ROIs' anatomical contours. Voxel based regression analysis A voxel-wise ordinary least squares (OLS) multiple regression analysis was performed using the AFNI (http://afni.nimh.nih.gov/afni) and SPM5 (Wellcome Department of Imaging Neuroscience, London, UK) software packages. Pre-processing included the following steps for each subject: the functional volumes were corrected for differences in slice acquisition times and head motion, and subsequently smoothed with an 8-mm FWHM Gaussian kernel. In addition to the gamma-convolved pain rating curve, representing the effect of interest, the linear model for the individual subjects' fMRI data analysis included the motion parameters obtained from the preprocessing motion-correction step and the global fMRI signal (average of the fMRI signal time courses over all in-brain voxels), as confounds. The map representing the estimated coefficient for the psychophysical pain rating regressor (the “beta” image) was scaled voxel-wise to represent percent change with respect to the voxel temporal mean. The beta image was then spatially warped to the standard Montreal Neurological Institute (MNI) space by applying the set of spatial transformation parameters estimated by applying the SPM5 segmentation routines to the high-resolution T1-weighted anatomical image, previously registered to the EPI dataset; the warped beta image was resampled to a 3 × 3 × 3 mm voxel size. The group analysis was performed as a one-sample t-test on the MNI-warped beta images, to provide a (random-effect) population-level inference about the brain regions whose activity profile was significantly predicted by the pain rating curve. A double statistical threshold (single-voxel p b 0.01; minimum cluster size = 145 voxels) was adopted to achieve an experiment-wise significance level (corrected for multiple comparisons) of p b 0.05, as computed by Monte Carlo simulation with the AFNI routine AlphaSim (http://afni.nimh.nih.gov/pub/dist/doc/manual/AlphaSim.pdf). To compare the regression results with those of the multi-way PLS analysis, we performed one-sample t-tests to assess which of the 44 ROIs had an average beta value significantly different from zero at the

S. Favilla et al. / NeuroImage 90 (2014) 153–162

group level. Significance values were corrected for multiple comparisons using a two-stage false discovery rate approach weighted for ROI size (wFDR) as described by Benjamini and Heller (2007). This procedure controls the wFDR, i.e. the expected proportion of erroneously rejected ROIs out of all tested ROIs, at p b 0.05. The procedure first estimates the number of non-significant ROIs (stage 1) and then uses this information to enhance the statistical power of the correction (stage 2). A two-way (2 hemispheres × 22 ROIs) repeated-measure analysis of variance (ANOVA) was also performed, to test for possible differences in ROI-averaged beta values between the left and right hemisphere. N-partial least squares analysis Preprocessing steps Functional EPI volumes were corrected for differences in slice acquisition times and head motion. A multiple regression analysis was used to regress out the variance explained by the nuisance regressors (i.e., the six motion parameters and the global fMRI signal) from the EPI signal while preserving the variance explained by the regressor of interest, the gamma-convolved pain rating curve. The cleaned EPI volumes were resampled to 2 × 2 × 2 mm voxel size and L2-normalized (i.e., the sum of squares of each voxel's time series was set to 1). The first principal singular vector (1st SVD) of each subject's EPI time series was computed within each of the ROIs. The sign of the SVD vector was chosen so that the average of the (Fisher z-transformed) correlation coefficients of the SVD vector with each voxel time series within the ROI mask was positive. The output SVD vector was L2-normalized. Multi-way PLS (N-PLS) PLS regression aims to find a relationship between a set of predictor (independent) data, X ∈ RI × J × K, and a set of responses (dependent data) Y ∈ RI × M × N by maximizing XY covariance. Multi-linear PLS represents a generalization on partial least squares (PLS) to multi-way datasets (Bro, 1996) applied to the study of an N-way problem in different application areas. Following the standard notation in the literature on multi-way statistical approaches (Kiers, 2000) we denote scalars and vectors by lower italics (a) and bold case (x) letters, respectively; bold-face capitals are used for matrices (X) and underlined bold-face capitals for three-way arrays (X). In general, the independent X and dependent Y variables are decomposed in a way such that the scores vectors from this model have pair-wise maximal covariance (de Jong, 1998; Smilde, 1997) through a regression model. The estimates of first mode scores matrix (T), second and third mode weights matrices (WJ, WK), and the regression coefficients (b) obtained by the N-PLS are generally satisfactory for the computation of Y but may be unsatisfactory in the fit of X. An improved version of the original N-PLS method (see detailed description in Bro et al. (2001) overcomes this drawback estimating a core array, G, based on the already estimated N-PLS scores and weights. In the improved N-PLS, the X and Y arrays of independent and dependent data, respectively, are modeled by the so-called Tucker3 decomposition (for further details, see Appendix A.a in Supplementary material), which enables the interaction of each component in each mode with all components in other modes. N-PLS regression is therefore a sequential algorithm where one component is computed at a time and each of these components is calculated by first finding the variable weights (W) of that component. By regressing the data onto their weight vectors, a score vector is found in the X-space providing a least squares model of the X data. Finally, by choosing the weights such that the covariance between X and Y is maximized, a regression model is obtained. In this study, we defined the independent variables block X ∈ RI × J × K as a three-way array of BOLD-fMRI signals (first principal singular vector (1st-SVD) of the BOLD signals belonging to a given ROI (representing subsequent time intervals) × n° ROIs × n° of volunteers)

155

sized 270 × 44 × 17 (Fig. 1), and the dependent variables block Y ∈ RI × K as the matrix (sized 270 × 17) containing the recorded psychophysical profile of perceived pain intensity for each subject. The choice of defining time as mode #1 was motivated by the aim of our study: namely, we wanted to identify those ROI time series which were strongly associated (i.e., in terms of covariance) with the psychophysical pain profile of each volunteer. We considered only the 270 signal time points common to all subjects, spanning the time period from 4 min before until 14 min after the ascorbic acid injection. Diagnostics and validation A generally accepted criterion to assess the number of components (latent variables) for a multi-way model is to select the dimensionality, i.e., the number of latent variables to retain, corresponding to the minimum root-mean-square error of cross-validation (RMSECV) (Smilde et al., 2004), which is a measure of the model ability to predict entries that were not used to build the model and is defined as follows: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u I  2 u1 X JK ^i −yJK RMSECV ¼ t y imeas I i¼1 ^JK where I is the number of time points, y i are the values of the predicted data points that were not included in the model formulation, and y JK i_meas are the values of the measured data points for each ROI and for each subject. Venetian blind cross-validation with 10 segmentations was used, thus in each cross-validation iteration 27 time points were left out and estimated by the model built on the remaining ones. The identification of individual samples and/or variables which are most influential for a given model is important and it may also provide useful information concerning the presence of outliers. Leverages have been mainly developed for linear regression problems (Hoaglin and Welsch, 1978; Smilde et al., 2004) and can be used in multilinear regression analysis for assessing the relative importance of different objects or variables according to the model. In general, a high value indicates an influential sample or variable, whereas a low value indicates the opposite (see Appendix A.b in Supplementary Material). In a regression model, the main goal is to make the residual Ey of Y small. Diagnostics based on Y-residual are well described in standard regression literature (see for instance, Martens and Næs, 1989). In order to identify which ROIs mainly contribute to the achievement of good estimation or mainly influence the given model, as well as to

Fig. 1. Three-way data set arrangement of X (I × J × K): Mode 1, time points (I = 270 TR); Mode 2, first principal singular vector (1st-SVD) of the BOLD-fMRI signals belonging to a region of interest (ROI) in subsequent time intervals (J = 44 ROIs); Mode 3, volunteers (K = 17).

156

S. Favilla et al. / NeuroImage 90 (2014) 153–162

evaluate the influence of individual volunteers, the residuals terms in Ex and Ey arrays were squared and summed across the different modes, thus providing the variation in the sum square of residuals (SSR) as a function of variables (ROIs) or volunteers. The behavior of the regression model was thus studied, for the second and third modes (corresponding to the ROIs and volunteers, respectively), in terms of leverage contribution (named hwj and hwk respectively) versus SSR. Small SSR and high leverage values indicate positive influence of ROIs/volunteers on the model, whereas higher SSR and small leverage values indicate negative influence on the model. The performance of the model was evaluated by measuring the rel  ^ −ymeas k2 =kymeas k where y ^ and ymeas ative reconstruction error ρ2 ¼ ky are the vectors of the calculated and real pain intensities of one subject, respectively. The results of the N-PLS model are presented in Supplementary Table 1. The root mean squares error in cross-validation (RMSECV) does not show a minimum but slowly decreases, and models with more than two latent variables (LVs) do not show a significant decrease in cross-validation error, according to an F-test at 95% confidence level. Here, N-PLS is mainly used as a projection method to obtain LVs co-varying with the psychophysical response, rather than to exactly reproduce it, since our main goal was to identify those ROI time series that are most involved in determining the subjective pain perception profile. Thus, we consequently decided to retain four LVs as a good compromise between a low number of latent variables and a satisfactory percentage of variance explained for Y. The N-PLS analysis was carried out using the freely available N-way 3.1 (Andersson and Bro, 2000).

Variable importance in projection: extended algorithm to three-way dataset Variable importance in projection (VIP) (Chong and Jun, 2005; Wold et al., 1993) is a score value useful to summarize the contribution of each variable to the model. VIP coefficients thus represent the importance of each x variable in accounting for both the x- and the y-variance. Since the y-variances are predicted from the x-variances, VIP allows ranking the X variables according to their explanatory power of Y. The VIP parameter is based on the idea of factor models. In linear factor models, several components summarizing the data are extracted either to explain the variance in the data, e.g., as in PCA, or to capture the correlation between two datasets, e.g., as in PLS or in canonical correlation analysis. The components extracted in these linear factor models are linear combinations of the data variables. For each variable x in the independent array X, a VIP value is computed quantifying the importance of x in explaining the response Y. The VIP value depends both on the correlation between X and Y, and on the relative contribution of each latent variable in explaining the variance of Y (see Appendix A.c in Supplementary materials, for mathematical description of VIP). The original code is suitable for two-way regression analysis (Wold et al., 1993). In this study, we applied an extension of the original two-way code of VIP to multi-way arrays (3w-VIP2) (Favilla et al., 2013), in order to highlight the independent variables (i.e., ROIs' BOLD

time courses) that are most relevant to describe the properties of interest (i.e., the time course of perceived pain). As in any ranking method, the choice of the VIP threshold is a critical issue. Since the average of squared VIP scores equals 1, we followed the general criterion to select those variables with squared VIP score greater than 1 (Wold et al., 1993): because they have an above average influence on the model, they are considered the most relevant for explaining Y.

Results Psychophysics In all subjects, the ascorbic acid injection induced a moderate to strong pain lasting several minutes. The time course of the perceived pain intensity was quite variable across subjects (Supplementary Fig. 2), confirming previous observations (Porro et al., 1998). After a common period of sharp and localized burning pain dominating in the first tens of seconds, the i.m. injection usually produced a dull or cramp-like pain sensation, whereas the s.c. injection was usually described as burning. The mean time course of the perceived pain intensity for subjects who received the i.m. ascorbic acid injection (n = 9) or the s.c. ascorbic acid injection (n = 8) is shown in Supplementary Fig. 3. No between-group significant differences were observed with regard to either the peak pain intensity (mean ± SD on a 0–100 scale: i.m. 48.5 ± 31.7; s.c. 51.1 ± 23.1; t(15) = −0.19; p = 0.849) or the area under the curve of pain intensity ratings across the whole experiment (mean ± SD: i.m.: 3775 ± 2875; s.c.: 4356 ± 1950; t(15) = −0.48; p = 0.638).

Regression analysis Results of the whole-brain multiple regression analysis are shown in Table 1. We found three clusters – with peaks located in the insular cortex (bilaterally) and in the anterior mid-cingulate cortex (Fig. 2) – whose activity profiles were significantly correlated with the individual pain rating curve. Results for the ROI-averaged regression beta values are shown in Table 2 and Fig. 3, top. Activity in 9 of the 44 selected ROIs displayed a significant relationship with the individual psychophysical profiles (namely, an average beta value significantly different from zero, corrected for multiple comparisons): the bilateral posterior insula (pINS), parietal operculum (PO), mid-cingulate cortex (MCC) and putamen (Put) showed positive beta values, whereas the ipsilateral (left hemisphere) precuneus showed negative beta values. The bilateral anterior insula (aINS) and the contralateral rostral medial superior frontal gyrus/paracentral lobule showed a trend towards positive beta values (uncorrected p b 0.05), but failed to reach corrected significance. ANOVA yielded no significant differences in ROI-averaged beta values between the two hemispheres (hemisphere × ROI interaction: F21,336 = 0.69; p = 0.839).

Table 1 Whole brain voxel-based regression: foci correlated with the individual pain rating profile. Anatomical region

Left insula Left insula Right insula Right insula Anterior mid-cingulate (bilateral)

BA

Cluster size (voxels)

Center of mass (mm) x

y

Peak (mm)

13

734

−39

−11

6

13

708

40

−6

10

24, 32

506

1

14

32

z

Peak p

x

y

z

(uncorrected)

−39 −39 42 33 −6 6

−24 −9 −3 −15 21 27

12 −3 9 15 21 21

0.000002 0.000003 0.000001 0.000002 0.000028 0.000132

Coordinates of the global center of mass and of the two most significant peaks are shown for each cluster. All coordinates refer to the MNI standard space. BA = Brodmann area.

S. Favilla et al. / NeuroImage 90 (2014) 153–162

157

Fig. 2. Voxel-wise regression group analysis: foci showing activity positively correlated with pain intensity, overlaid onto sagittal planes of the MNI template brain.

N-PLS analysis

Ranking of ROIs

The measured and fitted psychophysical profiles from three representative subjects are shown in Fig. 4. The low value of the averaged reconstruction error across all 17 subjects (ρ2 = 0.25 ± 0.11, mean ± SD) represents a good evidence of model consistency in reproducing different perceptual behaviors, starting from the individual fMRI data. Scores plot for the first mode (the temporal mode) of X N-PLS model are shown in Supplementary Fig. 4. The behavior of the regression model can be further studied in terms of leverage contribution for the second (hwj) and third (hwk) modes, corresponding to the ROIs and volunteers, respectively. Fig. 5 reports the scatter plot of leverage for mode #2 (ROIs) versus the residual sum of squares (SSR) for the X array. ROIs whose signals have been well modeled (small SSR) and that, at the same time, contribute to the model with unique information (high hwj) include the bilateral aINS and MCC, and the ipsilateral PO. Overall, leverage analysis for mode #2 (ROIs; Fig. 5) and #3 (volunteers; Supplementary Fig. 5) versus the residual sum of squares (SSR) did not reveal anomalous cases in term of either ROIs or subjects. Furthermore, a relatively random distribution of the Y residuals versus volunteers suggests a good fit of the model; no outliers were identified (Supplementary Fig. 6).

Due to the high efficiency of the model, we tested the capability of the 3w-VIP2 parameter to identify the brain areas most involved in the estimation of the psychophysical pain intensity profile. Twelve ROIs showed a 3w-VIP2 value N 1 (Table 3 and Fig. 3, bottom). The averaged BOLD signal changes of these 12 ROIs are shown in Fig. 6, in comparison with the mean psychophysical profile of the 17 volunteers. A good correspondence was found between BOLD time profiles of the most relevant ROIs and the pain intensity profile throughout the experiment. Notably, the bilateral mid-cingulate, posterior insular cortex and parietal operculum appeared to be more specifically involved during the ascending phase of pain. The BOLD time profiles of the 6 ROIs with the lowest 3w-VIP2 values did not show any changes related to the psychophysical profiles (Fig. 6, bottom). Eight of the ROIs with a 3w-VIP2 value N1, i.e. the bilateral pINS, PO, MCC and Put, also showed averaged beta values significantly greater than zero in the regression analyses (see Table 2 and Fig. 3). Discussion We applied multi-variate and multi-way analyses to the fMRI data acquired during ascorbic acid nociceptive stimulation in healthy

Table 2 Brain regions related to the individual pain rating profile (ROI-based regression analysis). Anatomical region

Hemisphere

Mean norm. beta

t

p (uncorr.)

p (wFDR)

Left posterior insula Right posterior insula Right parietal operculum Right mid-cingulate cortex Left parietal operculum Left mid-cingulate cortex Right putamen Left putamen Left anterior insulaa Right anterior insulaa Right medial superior frontal gyrus/paracentral lobulea Left superior parietal lobulea Left precuneus

I C C C I I C I I C C I I

0.272 0.226 0.432 0.289 0.304 0.281 0.167 0.194 0.292 0.220 0.125 −0.122 −0.306

7.021 4.632 3.774 3.742 3.507 3.491 3.206 2.904 2.683 2.613 2.469 −2.351 −2.835

0.00000 0.00028 0.00166 0.00178 0.00292 0.00302 0.00551 0.01035 0.01635 0.01885 0.02519 0.03186 0.01194

0.00008 0.00389 0.02073 0.01418 0.02157 0.01657 0.02680 0.04517 0.05339 0.05217 0.06554 0.07306 0.04717

Only ROIs showing mean normalized beta values different from 0 are shown. wFDR = cluster-size weighted false discovery rate. C and I = contralateral and ipsilateral to the injected hand, respectively. a Regions with p b 0.05 uncorrected, but not significant for wFDR.

158

S. Favilla et al. / NeuroImage 90 (2014) 153–162

volunteers. An advantage of this model is that it shares several features of clinical inflammatory pain, such as prolonged stimulation of deep structures and fluctuations of ongoing pain over a relatively wide range of the perceived intensity, in the absence of simultaneous exteroceptive input. Changes of perceived pain intensity usually occur over seconds, a time frame compatible with BOLD-fMRI (Bandettini and Ungerleider, 2001), as is the case in clinical pain (Baliki et al., 2006; Geha et al., 2007). Furthermore, while in experimental paradigms based on brief noxious stimuli abrupt changes in attention and arousal are entirely time-locked to pain perception, in our case the long duration of the pain sensation allows for a better decoupling of these processes. The N-PLS modeling combined with the novel three-way extension of the VIP method allowed us to identify a group of brain regions whose activity is more strongly related to the psychophysical pain intensity profile. To our knowledge, this approach has never been applied to the analysis of the human pain system, and more generally to fMRI data. Three main conclusions can be drawn. First, both sub-cortical and cortical structures appear to be involved in coding for subjective pain intensity. Second, most pain-related regions appear to be almost equally active in the two hemispheres. Third, the most influential cortical regions are the mid-cingulate and anterior insula, whereas no significant contribution of anterior parietal areas could be detected. A similar array of brain areas could be identified by ROI-based analysis of average regression beta-values. Methodological aspects

Fig. 3. Mean normalized beta values obtained from the voxel-wise regression analysis (top) and 3w-VIP2 scores (denoted as VIP in the figure) values (bottom), for regions showing a significant positive effect in at least one of the two analyses (mean beta value N 0 or VIP2 score N 1). C and I: contralateral and ipsilateral to the injected hand, respectively. aINS: anterior insula; pINS: posterior insula; MCC: mid-cingulate cortex; MSFG: medial superior frontal gyrus; MT: medial thalamus; PO: parietal operculum; Put: putamen.

ρ =0.17

Volunteer 10

1

psychophysical pain profile

Task In the present study, subjects performed a continuous task of explicit pain intensity coding, involving not only the allocation of attentional resources for monitoring nociceptive inflow, but also working memory and decision-making processes. However, the time profiles of these latter processes are likely to be relatively constant throughout the pain period, and therefore very different from those of pain (which featured a sharp rise at the onset of stimulation, followed by a slow decline in the majority of our subjects). Another potential confounding factor is the motor activity associated with turning the knob of the rating device. We do not believe this to be a significant contributor to our results, given the likely discrepancy in the respective time courses of pain intensity and motor activity (see Supplementary Fig. 7). Our interpretation is further supported by the results of GLM analysis, which did not

2

0.5 0 1

2

3

4

5

6

7

8

1

9

10

11

12

Volunteer 13

13

14

15

16

17

16

17

16

17

ρ =0.20 2

0.5 0 1

2

3

4

5

6

7

8

1

9

10

11

12

Volunteer 15

13

14

15

ρ =0.33 2

0.5 0 0

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

18

Time (min) Fig. 4. Comparison between measured psychophysical data (in red) and reconstructed pain intensity profiles from the N-PLS model (in green), in three representative volunteers. The corresponding reconstruction errors ρ2 are also shown.

S. Favilla et al. / NeuroImage 90 (2014) 153–162

159

26: MCC − I 0.2

Leverage2 (hwj)

25: MCC − C 18: PO − I 21: aINS − C 22: aINS − I

0.15

0.1

0.05

0

5

6

7

8

9

10

11

12

13

14

15

SSR2 Fig. 5. Leverage plotted against sum-squares of residuals (SSR) for the second mode (ROIs) of X N-PLS model. Values for each of the 44 ROIs are shown. The bilateral anterior insula (aINS) and mid-cingulate cortex (MCC), and the ipsilateral parietal operculum (PO) have large hwj and small residuals (SSR2), indicating that they positively contribute to the model.

show any significantly correlated activity with the pain rating profile within primary motor areas, even at a very liberal statistical threshold (Supplementary Fig. 8), and by the N-PLS results, where primary motor areas displayed very low VIP scores. In previous studies using the same model, both changes in voluntary motor output (ballistic upper limb extension) (Bonifazi et al., 2004) and vegetative reactions (Porro et al., 2003) showed a shorter time frame than the perceived pain. For instance, heart rate changes are evident immediately after the onset of pain following ascorbic acid (Porro et al., 2003) or hypertonic saline injection (Burton et al., 2009), but rapidly decline afterwards, usually returning to baseline within 60 s of the injection. On these grounds, we believe that the described pattern of brain activity is indeed related mainly to pain perception.

In the present study, N-PLS was used to model the psychophysical pain profile in different subjects on the basis of BOLD-fMRI signals. When compared with the ROI-based regression analysis, this approach presents some advantages. First, inter-subject variability is naturally accounted for in the model, without the need for a two-stage estimation process, as in the case of a GLM-based group analysis. Second, it is easy to recover information about characteristic psychophysical profiles in the subjects group, and the regions that are chiefly involved. Finally, a ranking of the most relevant ROIs is possible by using a parameter (the VIP score) that is inherent to the model and ensures that the emphasis is placed on sub-arrays of variables, instead of on single variables, taking into account the correlation among them. Brain system involved in pain intensity coding

Multi-way data analysis In recent years, multi-way models have been successfully applied in a wide variety of fMRI studies (Andersen and Rayens, 2004; Barnathan et al., 2011; Martinez-Montes et al., 2004; Mørup et al., 2008). Specifically, the advantage of N-PLS is that it simultaneously decomposes the fMRI data and the design or task matrix (or array) in a way that maximizes the covariance between them. This means that, whereas decomposition methods, such as principal components analysis or independent components analysis find components that account for the highest variance in the fMRI data alone, PLS finds components that best explain the relation between the data and the design (Krishnan et al., 2011).

Table 3 Brain regions related to the individual pain rating profile (N-PLS analysis). Anatomical region

Hemisphere

3w-VIP2 score

Left mid-cingulate cortex Right anterior insula Left anterior insula Right mid-cingulate cortex Left parietal operculum Left putamen Right posterior insula Left posterior insula Right parietal operculum Right medial superior frontal gyrus/paracentral lobule Right putamen Left medial thalamus

I C I C I I C I C C

3.709 3.695 3.558 3.547 3.165 3.092 2.488 2.431 1.918 1.730

C I

1.698 1.209

C and I = contralateral and ipsilateral to the injected hand, respectively. Only areas with a 3w-VIP2 score N1 are listed.

The array of cortical areas revealed by our present approaches closely matches the one which has been recently shown to be involved in the prediction of cutaneous thermal pain in healthy volunteers, using a machine-learning regression technique (Wager et al., 2013); specifically, the anterior and posterior insula and mid-cingulate cortex ranked high in our VIP analysis and significantly contributed to pain prediction in the study of Wager et al. (2013). The identification of a similar pattern of brain areas, across different experimental models (involving noxious stimulation of superficial or deep tissues) and data analysis approaches, strongly supports the involvement of these regions in pain perception. The posterior insula and the adjoining parietal operculum is the main sensory receiving area of the spinothalamic system in primates, and is considered a key region for pain generation in humans on the basis of functional and brain lesion studies (Garcia-Larrea, 2012; Oertel et al., 2012; Treede et al., 2000). The anterior insula has long been associated with pain (Peyron et al., 2000) and likely plays an important role for the assessment of stimulus intensity (Baliki et al., 2009; Wiech et al., 2010). More generally, it has been hypothesized to be the neural substrate for awareness of feelings from the body across time (Craig, 2009, 2011) and, together with the mid-cingulate cortex, has been shown to be strongly implicated in autonomic control (Critchley, 2005). It is interesting to note that, in our experiment, the parietal operculum exhibited the fastest reaction to the stimulus and thus appeared to be more specifically involved in the initial sharp ascending phase of the pain sensation. The aMCC and pINS time courses, on the other hand, seemed to be a better match overall for the entire span of the psychophysical profile. All of the above-mentioned cortical regions are commonly considered to be part of the pain system. Activity in the paracentral lobule, likely corresponding to the supplementary motor area, has been also associated with nociceptive processing (Duerden and Albanese, 2013;

160

S. Favilla et al. / NeuroImage 90 (2014) 153–162

Fig. 6. Comparison between the mean psychophysical profile and averaged BOLD signal changes in selected brain regions. From top to bottom: a) mean (±SD, n = 17) time course of the perceived pain intensity in all volunteers; b) averaged BOLD signal changes in the 6 brain regions showing the highest 3w-VIP2 scores (N3); c) averaged BOLD signal changes in the 6 brain regions showing 3w-VIP2 scores b3 but N1; d) averaged BOLD signal changes in the 6 brain regions showing the lowest 3w-VIP2 scores (b0.1). pACC: (perigenual) anterior cingulate cortex; LT: lateral thalamus; MPF: medial prefrontal cortex; rMSFG: rostral medial superior frontal gyrus; PostCG-Mid: post-central gyrus, middle. C and I: contralateral and ipsilateral to the injected hand, respectively. Other abbreviations as in Fig. 3.

Oshiro et al., 2009; Peyron et al., 2000), although its precise role awaits further demonstration. Overall, these findings confirm the hypothesis that cognitive representation of pain is attributable mainly to the activity of frontal, cingulate and insular regions (Fairhurst et al., 2012; Raij et al., 2009). The top-ranking brain regions identified here were cortical areas. However, the medial thalamus, and especially the putamen showed activity significantly related to pain intensity. Medial thalamic nuclei are known to encode stimulus intensity (Lenz et al., 2004) and to send nociceptive information to the mid- and anterior cingulate cortex; these regions are consistently activated by noxious stimuli (Duerden and Albanese, 2013), likely play a key and complex role in pain perception (Buchel et al., 2002; Oshiro et al., 2009; Vogt, 2005), and might be involved in top-down recurrent projection from the cortex to the thalamus. The basal ganglia have been shown to be active during prickle sensations evoked by noxious cold stimuli (Davis et al., 2002) and during intensity discrimination of noxious heat (Oshiro et al., 2009). Nociceptive input is known to be bilaterally distributed in the brain (see Olausson et al., 2001). A recent meta-analysis suggested a right hemisphere predominance in pain processing (Duerden and Albanese, 2013), especially in the insula and cingulate cortex, whereas we found no clear evidence for such hemispheric bias. The postcentral gyrus was not identified in any of our group analyses, whereas conflicting results have been obtained in other studies on the neural system predicting thermal pain (Brodersen et al., 2012; Wager et al., 2013). At the individual level, clusters related to pain intensity

coding in the primary somatosensory cortex (SI) have been found to be small and to show large interindividual variations in their extent and precise location (Porro et al., 1998, 2002). Previous event-related fMRI studies using short-term cutaneous noxious stimuli have also shown a predominant involvement of anterior mid-cingulate and insular/opercular regions during the intensity coding phase, whereas the anterior parietal cortex was mainly involved during the pain perception phase (Kong et al., 2006; Lui et al., 2008) and in short-term memory of pain (Albanese et al., 2007). Limitations of the study Both the N-PLS modeling and the averaged beta values analyses were based on a relatively limited number of pre-selected anatomical ROIs, chosen on strong prior assumptions about their potential involvement in somatosensory processing and pain perception. Gray matter ROIs were drawn in the native brain space in individual subjects, to account for between-subject variability in sulcal and gyral location; this approach has been shown to provide more accuracy to detect an experimental effect, relative to group or atlas-derived ROIs (Brett et al., 2002; Mitsis et al., 2008). However, the selected ROIs did not cover the whole brain, and therefore the contribution of other cortical areas which have been shown to be involved in pain processes, such as the lateral aspect of the frontal pole (Coghill et al., 2003; Derbyshire et al., 1997) and dorso-lateral prefrontal cortex (Oshiro et al., 2009), may have been missed.

S. Favilla et al. / NeuroImage 90 (2014) 153–162

Furthermore, the use of s.c. and i.m. injections in different subjects prevented us from detecting specific neural signatures associated with stimulation of different tissues, beyond their common pattern. Future directions Our present findings demonstrate that, across individual subjects, a complex array of brain regions show activity changes over time that are related to the perceived changes in pain intensity. More complex N-PLS models on brain clusters derived from data-driven parcellation methods (Craddock et al., 2012) can be foreseen. Future studies using these multiway approaches on fMRI data obtained in different experimental models of pain, and in pain patients, will hopefully help to understand the neural bases of pain perception at the individual level (Coghill, 2010; Tracey, 2011), thus paving the way to specific and more powerful therapeutic approaches. Acknowledgments The authors thank the Fondazione Cassa di Risparmio di Modena (FCRM) for its financial support to the Modena MR center. This work was funded by a grant from FCRM — International Projects to C.A.P. Appendix A. Supplementary data Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.neuroimage.2014.01.001. References Albanese, M.C., Duerden, E.G., Rainville, P., Duncan, G.H., 2007. Memory traces of pain in human cortex. J. Neurosci. 27, 4612–4620. Andersen, A.H., Rayens, W.S., 2004. Structure-seeking multilinear methods for the analysis of fMRI data. Neuroimage 22, 728–739. Andersson, C.A., Bro, R., 2000. The N-way Toolbox for MATLAB. Chemom. Intell. Lab. Syst. 52, 1–4. Apkarian, A.V., Bushnell, M.C., Treede, R.D., Zubieta, J.K., 2005. Human brain mechanisms of pain perception and regulation in health and disease. Eur. J. Pain 9, 463–484. Baliki, M.N., Chialvo, D.R., Geha, P.Y., Levy, R.M., Harden, R.N., Parrish, T.B., Apkarian, A.V., 2006. Chronic Pain and the Emotional Brain: Specific Brain Activity Associated with Spontaneous Fluctuations of Intensity of Chronic Back Pain. J. Neurosci. 26, 12165–12173. Baliki, M.N., Geha, P.Y., Apkarian, A.V., 2009. Parsing pain perception between nociceptive representation and magnitude estimation. J. Neurophysiol. 101, 875–887. Bandettini, P.A., Ungerleider, L.G., 2001. From neuron to BOLD: new connections. Nat. Neurosci. 4 (9), 864–866. Barnathan, M., Megalooikonomou, V., Faloutsos, C., Faro, S., Mohamed, F.B., 2011. TWave: high-order analysis of functional MRI. Neuroimage 58, 537–548. Benjamini, Y., Heller, R., 2007. False discovery rates for spatial signals. J. Am. Stat. Assoc. 102, 1272–1281. Bonifazi, M., Ciccarone, G., della, V.R., Spidalieri, R., Rossi, A., 2004. Influences of chemically-induced muscle pain on power output of ballistic upper limb movements. Clin. Neurophysiol. 115, 1779–1785. Brett, M., Johnsrude, I.S., Owen, A.M., 2002. The problem of functional localization in the human brain. Nat. Rev. Neurosci. 3, 243–249. Bro, R., 1996. Multiway calibration. Multilinear PLS. J. Chemom. 10, 47–61. Bro, R., Smilde, A.K., de Jong, S., 2001. On the difference between low-rank and subspace approximation: improved model for multi-linear PLS regression. Chemom. Intell. Lab. Syst. 58, 3–13. Brodersen, K.H., Wiech, K., Lomakina, E.I., Lin, C.S., Buhmann, J.M., Bingel, U., Ploner, M., Stephan, K.E., Tracey, I., 2012. Decoding the perception of pain from fMRI using multivariate pattern analysis. Neuroimage 63, 1162–1170. Buchel, C., Bornhovd, K., Quante, M., Glauche, V., Bromm, B., Weiller, C., 2002. Dissociable neural responses related to pain intensity, stimulus intensity, and stimulus awareness within the anterior cingulate cortex: a parametric single-trial laser functional magnetic resonance imaging study. J. Neurosci. 22, 970–976. Burton, A.R., Birznieks, I., Bolton, P.S., Henderson, L.A., Macefield, V.G., 2009. Effects of deep and superficial experimentally induced acute pain on muscle sympathetic nerve activity in human subjects. J. Physiol. 587, 183–193. Chong, I.G., Jun, C.H., 2005. Performance of some variable selection methods when multicollinearity is present. Chemom. Intell. Lab. Syst. 78, 103–112. Coghill, R.C., 2010. Individual differences in the subjective experience of pain: new insights into mechanisms and models. Headache 50, 1531–1535. Coghill, R.C., McHaffie, J.G., Yen, Y.F., 2003. Neural correlates of interindividual differences in the subjective experience of pain. Proc. Natl. Acad. Sci. U. S. A. 100, 8538–8542.

161

Craddock, R.C., James, G.A., Holtzheimer III, P.E., Hu, X.P., Mayberg, H.S., 2012. A whole brain fMRI atlas generated via spatially constrained spectral clustering. Hum. Brain Mapp. 33, 1914–1928. Craig, A.D., 2009. Emotional moments across time: a possible neural basis for time perception in the anterior insula. Philos. Trans R. Soc. Lond. B Biol. Sci. 364, 1933–1942. Craig, A.D., 2011. Significance of the insula for the evolution of human awareness of feelings from the body. Ann. N. Y. Acad. Sci. 1225, 72–82. Critchley, H.D., 2005. Neural mechanisms of autonomic, affective, and cognitive integration. J. Comp. Neurol. 493, 154–166. Davis, K.D., Pope, G.E., Crawley, A.P., Mikulis, D.J., 2002. Neural correlates of prickle sensation: a percept-related fMRI study. Nat. Neurosci. 5, 1121–1122. de Jong, S., 1998. Regression coefficients in multilinear PLS. J. Chemom. 12, 77–81. Derbyshire, S.W., Jones, A.K., Gyulai, F., Clark, S., Townsend, D., Firestone, L.L., 1997. Pain processing during three levels of noxious stimulation produces differential patterns of central activity. Pain 73, 431–445. Duerden, E.G., Albanese, M.C., 2013. Localization of pain-related brain activation: a metaanalysis of neuroimaging data. Hum. Brain Mapp. 34, 109–149. Fairhurst, M., Fairhurst, K., Berna, C., Tracey, I., 2012. An fMRI study exploring the overlap and differences between neural representations of physical and recalled pain. PLoS One 7, e48711. Favilla, S., Durante, C., LI Vigni, M., Cocchi, M., 2013. Assessing features relevance in NPLS models by VIP. Chemom. Intell. Lab. Syst. http://dx.doi.org/10.1016/j.chemolab.2013. 05.013. Garcia-Larrea, L., 2012. The posterior insular–opercular region and the search of a primary cortex for pain. Neurophysiol. Clin. 42, 299–313. Geha, P.Y., Baliki, M.N., Chialvo, D.R., Harden, R.N., Paice, J.A., Apkarian, A.V., 2007. Brain activity for spontaneous pain of postherpetic neuralgia and its modulation by lidocaine patch therapy. Pain 128, 88–100. Hoaglin, D., Welsch, E.R., 1978. The hat matrix in regression and ANOVA. Am. Stat. 32, 17–22. Kiers, H.A.L., 2000. Towards a standardized notation and terminology in multiway analysis. J. Chemom. 14, 105–122. Kong, J., White, N.S., Kwong, K.K., Vangel, M.G., Rosman, I.S., Gracely, R.H., Gollub, R.L., 2006. Using fMRI to dissociate sensory encoding from cognitive evaluation of heat pain intensity. Hum. Brain Mapp. 27, 715–721. Krishnan, A., Williams, L.J., McIntosh, A.R., Abdi, H., 2011. Partial least squares (PLS) methods for neuroimaging: a tutorial and review. Neuroimage 56, 455–475. Kulkarni, B., Bentley, D.E., Elliott, R., Youell, P., Watson, A., Derbyshire, S.W., Frackowiak, R.S., Friston, K.J., Jones, A.K., 2005. Attention to pain localization and unpleasantness discriminates the functions of the medial and lateral pain systems. Eur. J. Neurosci. 21, 3133–3142. Legrain, V., Iannetti, G.D., Plaghki, L., Mouraux, A., 2011. The pain matrix reloaded: a salience detection system for the body. Prog. Neurobiol. 93, 111–124. Lenz, F.A., Ohara, S., Gracely, R.H., Dougherty, P.M., Patel, S.H., 2004. Pain encoding in the human forebrain: binary and analog exteroceptive channels. J. Neurosci. 24, 6540–6544. Lui, F., Duzzi, D., Corradini, M., Serafini, M., Baraldi, P., Porro, C.A., 2008. Touch or pain? Spatio-temporal patterns of cortical fMRI activity following brief mechanical stimuli. Pain 138, 362–374. Marquand, A., Howard, M., Brammer, M., Chu, C., Coen, S., Mourao-Miranda, J., 2010. Quantitative prediction of subjective pain intensity from whole-brain fMRI data using Gaussian processes. Neuroimage 49, 2178–2189. Martens, H., Næs, T., 1989. Multivariate Calibration. John Wiley and Sons, Chichester and New York. Martinez-Montes, E., Valdes-Sosa, P.A., Miwakeichi, F., Goldman, R.I., Cohen, M.S., 2004. Concurrent EEG/fMRI analysis by multiway partial least squares. Neuroimage 22, 1023–1034. McIntosh, A.R., Misic, B., 2013. Multivariate statistical analyses for neuroimaging data. Annu. Rev. Psychol. 64, 499–525. McIntosh, A.R., Bookstein, F.L., Haxby, J.V., Grady, C.L., 1996. Spatial pattern analysis of functional brain images using partial least squares. Neuroimage 3, 143–157. Mitsis, G.D., Iannetti, G.D., Smart, T.S., Tracey, I., Wise, R.G., 2008. Regions of interest analysis in pharmacological fMRI: how do the definition criteria influence the inferred result? Neuroimage 40, 121–132. Mørup, M., Hansen, L.K., Arnfred, S.M., Lim, L.H., Madsen, K.H., 2008. Shift-invariant multilinear decomposition of neuroimaging data. Neuroimage 42, 1439–1450. Oertel, B.G., Preibisch, C., Martin, T., Walter, C., Gamer, M., Deichmann, R., Lotsch, J., 2012. Separating brain processing of pain from that of stimulus intensity. Hum. Brain Mapp. 33, 883–894. Ohara, S., Crone, N.E., Weiss, N., Kim, J.H., Lenz, F.A., 2008. Analysis of synchrony demonstrates that the presence of “pain networks” prior to a noxious stimulus can enable the perception of pain in response to that stimulus. Exp. Brain Res. 185, 353–358. Olausson, H., Ha, B., Duncan, G.H., Morin, C., Ptito, A., Ptito, M., Marchand, S., Bushnell, M.C., 2001. Cortical activation by tactile and painful stimuli in hemispherectomized patients. Brain 124, 916–927. Oshiro, Y., Quevedo, A.S., McHaffie, J.G., Kraft, R.A., Coghill, R.C., 2009. Brain mechanisms supporting discrimination of sensory features of pain: a new model. J. Neurosci. 29, 14924–14931. Pardo, R., Helena, B.A., Cazurro, C., Guerra, C., Deban, L., Guerra, C.M., Vega, M., 2004. Application of two- and three-way principal component analysis to the interpretation of chemical fractionation results obtained by the use of the B.C.R. procedure. Anal. Chim. Acta. 523, 125–132. Paxinos, G., Mai, J.K., 2004. The Human Nervous System. Elsevier, San Diego. Peyron, R., Laurent, B., Garcia-Larrea, L., 2000. Functional imaging of brain responses to pain. A review and meta-analysis. Neurophysiol. Clin. 30, 263–288. Piche, M., Arsenault, M., Rainville, P., 2010. Dissection of perceptual, motor and autonomic components of brain activity evoked by noxious stimulation. Pain 149, 453–462.

162

S. Favilla et al. / NeuroImage 90 (2014) 153–162

Ploner, M., Lee, M.C., Wiech, K., Bingel, U., Tracey, I., 2010. Prestimulus functional connectivity determines pain perception in humans. Proc. Natl. Acad. Sci. U. S. A. 107, 355–360. Porro, C.A., 2003. Functional imaging and pain: behavior, perception, and modulation. Neuroscientist 9, 354–369. Porro, C.A., Cettolo, V., Francescato, M.P., Baraldi, P., 1998. Temporal and intensity coding of pain in human cortex. J. Neurophysiol. 80, 3312–3320. Porro, C.A., Baraldi, P., Pagnoni, G., Serafini, M., Facchin, P., Maieron, M., Nichelli, P., 2002. Does anticipation of pain affect cortical nociceptive systems? J. Neurosci. 22, 3206–3214. Porro, C.A., Cettolo, V., Francescato, M.P., Baraldi, P., 2003. Functional activity mapping of the mesial hemispheric wall during anticipation of pain. Neuroimage 19, 1738–1747. Prato, M., Favilla, S., Zanni, L., Porro, C.A., Baraldi, P., 2011. A regularization algorithm for decoding perceptual temporal profiles from fMRI data. Neuroimage 56, 258–267. Pravdova, V., Boucon, C., de Jong, S., Walczak, B., Massart, D.L., 2002. Three-way principal component analysis applied to food analysis: an example. Anal. Chim. Acta. 462, 133–148. Raij, T.T., Numminen, J., Narvanen, S., Hiltunen, J., Hari, R., 2009. Strength of prefrontal activation predicts intensity of suggestion-induced pain. Hum. Brain Mapp. 30, 2890–2897. Rossi, A., Decchi, B., 1997. Changes in Ib heteronymous inhibition to soleus motoneurones during cutaneous and muscle nociceptive stimulation in humans. Brain Res. 774, 55–61. Schulz, E., Zherdin, A., Tiemann, L., Plant, C., Ploner, M., 2012. Decoding an individual's sensitivity to pain from the multivariate analysis of EEG data. Cereb. Cortex 22, 1118–1123.

Smilde, A.K., 1997. Comments on multilinear PLS. J. Chemom. 11, 367–377. Smilde, A.K., Bro, R., Geladi, P., 2004. Multi-way Analysis with Applications in the Chemical Sciences. John Wiley & Sons, Ltd. Tracey, I., 2008. Imaging pain. Br. J. Anaesth. 101, 32–39. Tracey, I., 2011. Can neuroimaging studies identify pain endophenotypes in humans? Nat. Rev. Neurol. 7, 173–181. Treede, R.D., Apkarian, A.V., Bromm, B., Greenspan, J.D., Lenz, F.A., 2000. Cortical representation of pain: functional characterization of nociceptive areas near the lateral sulcus. Pain 87, 113–119. Vogt, B.A., 2005. Pain and emotion interactions in subregions of the cingulate gyrus. Nat. Rev. Neurosci. 6, 533–544. Wager, T.D., Atlas, L.Y., Lindquist, M.A., Roy, M., Woo, C.W., Kross, E., 2013. An fMRI-based neurologic signature of physical pain. N. Engl. J. Med. 368, 1388–1397. Wiech, K., Ploner, M., Tracey, I., 2008. Neurocognitive aspects of pain perception. Trends Cogn. Sci. 12, 306–313. Wiech, K., Lin, C.S., Brodersen, K.H., Bingel, U., Ploner, M., Tracey, I., 2010. Anterior insula integrates information about salience into perceptual decisions about pain. J. Neurosci. 30, 16324–16331. Wold, S., Johansson, E., Cocchi, M., 1993. 3DQSAR in Drug Design. In: Kubinyi, H. (Ed.). ESCOM Science Publishers, Leiden, Holland. Yousry, T.A., Schmid, U.D., Alkadhi, H., Schmidt, D., Peraud, A., Buettner, A., Winkler, P., 1997. Localization of the motor hand area to a knob on the precentral gyrus. A new landmark. Brain 120 (Pt 1), 141–157.