To date, relatively little is known about the spatiotemporal aspects of whole-brain blood oxygenation level-dependent (BOLD) responses to brief nociceptive stimuli. It is known that the majority of brain areas show a stimulus-locked response, whereas only some are characterized by a canonical hemodynamic response function. Here, we investigated the time course of brain activations in response to mechanical pain stimulation applied to participants' hands while they were undergoing functional magnetic resonance imaging (fMRI) scanning. To avoid any assumption about the shape of BOLD response, we used an unsupervised data-driven method to group voxels sharing a time course similar to the BOLD response to the stimulus and found that whole-brain BOLD responses to painful mechanical stimuli elicit massive activation of stimulus-locked brain areas. This pattern of activations can be segregated into 5 clusters, each with a typical temporal profile. In conclusion, we show that an extensive activity of multiple networks is engaged at different time latencies after presentation of a noxious stimulus. These findings aim to motivate research on a controversial topic, such as the temporal profile of BOLD responses, the variability of these response profiles, and the interaction between the stimulus-related BOLD response and ongoing fluctuations in large-scale brain networks.
Harmful stimuli applied to the skin activate small, myelinated Aδ-fibers and unmyelinated C-fibers. While the sensation mediated by Aδ-fibers is referred to as “first pain” and described as a sharp, pin prick, that mediated by C-fibers, referred to as “second pain”, is described as a dull, long-lasting, burning feeling. Previous functional magnetic resonance imaging (fMRI) studies have shown that the cortical response to painful contact heat is biphasic. Indeed, this biphasic response has been related to different possible physiological mechanisms such as first or second pain (Chen et al. 2002), or innocuous versus painful heat (Becerra et al. 2001; Wager et al. 2004; Moulton et al. 2005), but also nonpain-related processes such as participants' assessment of pain during stimulation (Moulton et al. 2005; Oshiro et al. 2007).
It is well established that painful stimuli trigger responses in a wide array of cortical and subcortical areas, including the primary and secondary somatosensory cortices (SI and SII), the cingulate and insular cortices, the thalamus, the prefrontal cortices, and the cerebellum (Peyron et al. 2000; Garcia-Larrea et al. 2003; Apkarian et al. 2005; Friebel et al. 2011); however, less research has been conducted on the temporal profile of the blood oxygenation level-dependent (BOLD) response to painful stimuli. A number of previous studies investigating the temporal profile of the BOLD response have focused either on the relationship between pain perception and the BOLD signal (Porro et al. 1998), or on the differences in the temporal properties of the hemodynamic signal in response to innocuous versus noxious stimuli in prespecified areas (Chen et al. 2002; Moulton et al. 2005). Fewer authors have proposed model-free or modified general linear model (GLM) whole-brain analysis of the temporal profile of the BOLD signal (Upadhyay et al. 2010; Moulton et al. 2012). For instance, Moulton and colleagues (Moulton et al. 2012) subdivided the BOLD response to heat stimuli into 3 temporal segments using time-shifted predictors in a GLM analysis. The results of this study indicated that only the late response correlates with perceptual ratings of intensity of the stimulus. Interestingly, the authors suggested a temporal differentiation between the brain responses to pain that likely reflects the involvement of different brain networks at different time latencies. Notably, most of the previous studies that have attempted to characterize the temporal profile of the BOLD response following painful stimuli have used stimuli of long duration (Chen et al. 2002; Moulton et al. 2005, 2012; Upadhyay et al. 2010). However, as evidenced by the finite impulse response theory (Oppenheim et al. 1983), only the use of brief stimuli is suitable for the effective characterization of the temporal profile of a response signal (Menon 2012). Moreover, in the majority of these studies, the possibility of investigating the whole-brain network response to noxious stimuli was spatially limited and not truly data driven.
The method most commonly employed to analyze fMRI data involves the creation of a set of predictors that are convolved with a prespecified BOLD response model (often a 2-gamma function; Glover 1999). This hypothesis assumes that the BOLD response is constant in time and space (Buckner et al. 1996). This method suffers from 2 principal problems (Gonzalez-Castillo et al. 2012): (1) A high number of false-negative findings in the presence of a low signal-to-noise ratio (Huettel and McCarthy 2001; Saad et al. 2003); (2) only areas showing the prespecified response are reliably detected. It has, in fact, recently been demonstrated that, under optimal conditions, almost all of the brain cortex (96%) shows a significantly stimulus-locked BOLD response (Saad et al. 2003; Gonzalez-Castillo et al. 2012). The majority of these areas (68%) show a response that is very different to the canonical hemodynamic response function (HRF). Other studies have emphasized that the BOLD response can be characterized by spatially variable delays and envelopes (Lai et al. 1993; Lee et al. 1995; Schacter et al. 1997; Buckner et al. 1998; Robson et al. 1998; Kruggel and von Cramon 1999a,b; Bandettini 2000; Saad et al. 2003; Handwerker et al. 2004; Thomason et al. 2005; Meltzer et al. 2008). Such variations have been attributed alternatively to hemodynamic delays, neuronal activity, anatomical differences, or fMRI noise (Weisskoff et al. 1993; Biswal et al. 1995; Biswal et al. 1996; Mitra et al. 1997; Kruggel and von Cramon 1999a,b; Saad et al. 2003; Chang et al. 2008).
Some promising alternatives for analyzing fMRI data without imposing any prespecified response model are represented by finite impulse response models (for more details see Goutte et al. 2000; Miezin et al. 2000; Ollinger et al. 2001), basis function (see for example Hossein-Zadeh et al. 2003; Woolrich et al. 2004), and temporal clustering techniques (Smolders et al. 2007).
In light of these previous suggestions proposing that painful stimuli are processed by different brain networks, each with a characteristic BOLD response, we aimed to characterize whole-brain responses to mechanical painful stimuli in a true data-driven fashion, that is, without imposing a predefined BOLD response or any a priori selection of regions of interest (ROIs). This methodology was employed to deal with the shortcomings of the traditional predictor-related fMRI analysis method (Glover 1999) under the hypothesis of a massive involvement in stimulus-locked brain areas during an active task (Gonzalez-Castillo et al. 2012). We investigated whether, following painful mechanical stimuli, different groups of areas (e.g. networks) activate (or deactivate) with a different temporal profile. We examined the time course of voxels of the whole-brain BOLD signal (Ploran et al. 2007) in response to mechanical painful stimuli (Magerl et al. 2001; Baumgartner et al. 2010) and used a hierarchical cluster analysis (Smolders et al. 2007; Sack et al. 2008; Gonzalez-Castillo et al. 2012) to classify the time course profile of voxels into clusters sharing similar temporal profiles. We assumed that employing data analysis that does not impose any predefined BOLD shape would enable us to demonstrate that the group of stimulus-locked areas that collaborate to process the experimental stimulus form a much wider picture than shown by previous studies. Subsequently, we specified the temporal characteristics of the signal in the detected areas using the BOLD latency mapping (BLM; Formisano and Goebel 2003) technique and event-related averages. BLM uses a linear fit method to characterize the parameters of the BOLD response as onset, duration, time to peak, and percentage of signal change (PSC). To establish the specificity of our findings, we also performed a control experiment in which nonpainful tactile stimuli were used.
Materials and Methods
Painful Stimulation Experiment
Seventeen healthy right-handed volunteers (8 females, mean age 28 ± 4.2 years) free of neurological or psychiatric disorders, not taking medications known to alter brain activity, and with no history of drug or alcohol abuse, participated in the study. Handedness was ascertained with the Edinburgh Inventory (Oldfield 1971). We obtained the written informed consent of each subject, in accordance with the Declaration of Helsinki. The study was approved by the institutional committee on ethical use of human subjects at the University of Turin.
Task and Image Acquisition
Prior to the fMRI session, the participants were instructed on the task and familiarized with the stimuli. The intensity of the stimuli was adjusted so as to obtain a painful stimulation in all subjects. We chose, for each participant, the stimulus that elicited the first painful sensation, using the modified method of limits (Greenspan and McGillis 1994). During the 4 runs of the experimental session, the participants were required to keep their eyes shut and relax while receiving the stimuli on their left and right hands. In this slow event-related paradigm, each run was composed of a sequence of 12 brief stimuli. The interstimulus interval was jittered pseudorandomly between 18 and 22 s to permit the complete envelope of the BOLD response and an adequate sampling of the hemodynamic function. No more than 3 stimuli were consecutively applied to the same hand. The location of the stimuli was changed slightly after each stimulus. The participants were required to use their fingers to indicate their rating of the mean perceived intensity of the stimuli at the end of the block. This procedure was adopted to avoid any spurious motor or cognitive response during the block. The stimuli were rated on a scale ranging from 0 (no pricking pain) to 10 (most intense pricking pain imaginable). The stimuli were applied using a hand-held 256-mN pin-prick probe with a flat cylindrical tip (diameter 250 µm). Mechanical pin-prick pain is primarily mediated by high-threshold mechanoreceptors and type I mechano-heat receptors (Magerl et al. 2001). Pin-prick stimulation has been used to evaluate mechanical pain in healthy (Baumgartner et al. 2010) and clinical (Rolke et al. 2006) populations. Previous fMRI studies have shown that areas activated in response to such mechanical painful stimuli are comparable with those following laser stimulations (Baumgartner et al. 2010). In addition, the brief duration of the stimulus (1 s) is particularly suited for temporal analysis of the BOLD signal, owing to the favorable ratio between the duration of the stimulus and the duration of the BOLD response to it (Oppenheim et al. 1983). Stimuli were applied by a trained experimenter; in this way the correct pressure and duration of the stimulus were guaranteed. The duration of the stimulus was ascertained using a metronome.
Images were acquired on a 1.5-T INTERA™ scanner (Philips Medical Systems) with a SENSE high-field, high-resolution head coil optimized for functional imaging. Functional T2*-weighted images were acquired using echoplanar imaging sequences, with a repetition time (TR) of 2000 ms, an echo time (TE) of 50 ms, and a 90° flip angle. The acquisition matrix was 64 × 64, with a 200-mm field of view (FoV). A total of 240 volumes were acquired, with each volume consisting of 19 axial slices, parallel to the anterior–posterior commissure (AC–PC); slice thickness was 4.5 mm with a 0.5-mm gap. In-plane resolution was 3.1 mm. Two scans were added at the beginning of functional scanning to reach a steady-state magnetization before acquiring the experimental data. The data from these scans were then discarded. Within a single session, a set of 3-dimensional (3D) high-resolution T1-weighted structural images was acquired for each participant, using a fast field echo sequence, with a 25-ms TR, an ultra-short TE, and a 30° flip angle. The acquisition matrix was 256 × 256, and the FoV was 256 mm. The set consisted of 160 contiguous sagittal images covering the whole brain. In-plane resolution was 1 mm × 1 mm and slice thickness 1 mm (1 × 1 × 1 mm3 voxels).
BOLD imaging data were preprocessed and analyzed using the BrainVoyager QX software (Brain Innovation, Maastricht, The Netherlands). Functional images were preprocessed as follows to reduce artifacts (Miezin et al. 2000): (1) The global intensity of the repeatedly measured images of a slice was corrected using a mean intensity adjustment. For each slice, we computed the average intensity across the first image, for each subsequent scan of the same slice, we computed the mean intensity and then scaled to obtain the same average slice intensity; (2) we corrected small head movements using 3D motion correction. All volumes were aligned spatially to the first volume by rigid body transformations, using a trilinear interpolation algorithm; (3) we employed a slice scan time correction to allow a whole volume to be treated as a single data point. The sequentially scanned slices comprising each volume were interpolated in time, using a signal sinc-interpolation algorithm; (4) we performed spatial data smoothing using a 3D Gaussian kernel with full-width at half-maximum (FWHM) of 8 mm; and (5) applied temporal filters to remove drifts due to scanner and physiological noise: Linear and nonlinear trends removal through a temporal high-pass filter eliminating frequencies <3 cycles in the time course were performed. For all temporal analyses, no temporal smoothing or low-pass filter was applied to preserve the temporal characteristics of the signal. Only for the GLM analysis, a temporal smoothing of 2.8-s FWHM was applied. After preprocessing, we followed a series of steps to obtain a precise anatomical location of brain activity to facilitate intersubject averaging. Each subject's slice-based functional scan was coregistered with their 3D high-resolution structural scan. The 3D structural dataset of each subject was then transformed into Talairach space (Talairach and Tournoux 1988): The cerebrum was translated and rotated into the AC–PC plane and the borders of the cerebrum were identified. Using the anatomical–functional coregistration matrix and the Talairach reference points, the functional time course of each subject was then transformed into Talairach space and the volume time course created.
To determine whether different groups of voxels had a specific time course in relation to the presentation of the stimulus, we applied fuzzy c-mean clustering (Smolders et al. 2007) to each voxel's time course in a time window of 22 s after the stimulus presentation. Fuzzy clustering has the advantage of performing a completely data-driven and unbiased parcellation of the voxels that exhibit a time course correlated with the stimulus presentation, that is, voxels showing a similar time course after the presentation of the stimulus are clustered together. This method decomposes the original fMRI time series into a predefined number of spatio-temporal modes, which include a spatial map and an associated cluster centroid time course. Group cluster maps were obtained using a 1-sample t-test; maps were computed at a statistical threshold of P < 0.05 corrected for multiple comparisons using false discovery rate correction (Benjamini and Hochberg 1995). Furthermore, the intersubject correspondence (ISC) index was calculated to show the exact percentage of subjects in which each cluster was found.
BOLD Latency Mapping
We estimated a series of parameters for single trials of the unsmoothed time courses (Richter et al. 2000; Formisano et al. 2002; Formisano and Goebel 2003) using the BLM plugin implemented in Brain Voyager (Esposito, personal communication). These parameters were: onset latency, time to peak, FWHM duration, latency and PSC (Lindquist and Wager 2007; Lindquist et al. 2009). With the use of the BLM technique, BOLD responses were estimated according to a piece-wise linear (trapezoidal) model fit of the slow event-related response (Richter et al. 2000; Formisano et al. 2002; Supplementary Fig. 1).
Confidence intervals were based on the variance of single-trial parameter estimates and the assumption of a t-distribution. This analysis can be performed either in a ROI-wise or in a voxel-wise fashion. The advantage of ROI-wise analysis is that it considers clusters as single entities, therefore allowing the statistical comparison of temporal features between clusters. In contrast, voxel-wise analysis enables the temporal parameters of single voxels within the clusters to be observed.
Specificity of the Clusters
An important question of this study was to ascertain whether, by decomposing the BOLD signal into clusters of voxels having a similar temporal profile, it would be possible to identify pain-specific activations. We used 2 approaches to answer this question. First, we combined meta-analytic tools and Bayesian inference techniques. We employed the Neurosynth Database for this purpose (Yarkoni et al. 2011). Neurosynth is a highly automated brain mapping database and framework that, using text mining and meta-analytic procedures, can be used to explore the representational characteristics of several neural and cognitive states. This framework can be used to explore the specificity of an observed pattern of activation given a search term. The problem of specificity calls for the reverse inference solution (Poldrack 2006). Indeed, the majority of neuroimaging studies provide a weak basis for determining what cognitive states a given brain pattern implies. Using Neurosynth, we quantified the forward inference or the probability that there would be activation in specific brain regions given the presence of a particular term P(activation|term), and the reverse inference or the probability that a term would occur in an article given the presence of activation in a particular brain region P(term|activation). By comparing these 2 analyses, we were able to create a map where commonly observed pain-related areas are coded with regard to their specificity for pain (Yarkoni et al. 2011). As a second step, we performed a control experiment in which innocuous brushing was used.
Nonpainful Stimulation (Control Experiment)
A second cohort of 12 healthy right-handed volunteers (6 females, mean age 29 ± 7.41 years) took part in the control experiment. One of the volunteers had previously taken part in the main experiment. The participants were free of neurological or psychiatric disorders, not taking medications known to alter brain activity, and had no history of drug or alcohol abuse. All were right-handed according to the Edinburgh Inventory (Oldfield 1971). We obtained the written informed consent of each participant, and the study was approved by the institutional committee on ethical use of human subjects at the University of Turin.
Task and Image Acquisition
In this control experiment, we used exactly the same design as before, except for the kind of stimulus. Here, we applied innocuous nonpainful tactile stimuli using a SenseLab™ Brush-05.
The stimuli were rated on a scale ranging from 0 (no brushing sensation) to 10 (most intense sensation imaginable). The stimulus was of the same duration as the pin-prick stimulation. As in the first experiment, the stimuli were applied to the dorsum of the hand. Images were acquired on a 1.5-T INTERA™ scanner (Philips Medical Systems) with a SENSE high-field, high-resolution head coil optimized for functional imaging. Functional T2*-weighted images and 3D high-resolution T1-weighted structural images were acquired using the same parameters used in the main experiment.
BOLD imaging data were preprocessed using the same preprocessing strategies as in the main experiment. After preprocessing, a series of event-related averages were created for each positive cluster detected with the temporal clustering strategy (see previous Materials and Methods section for details).
Perception of the Stimuli
The average rating of intensity for painful stimuli was 3.15 (±2); for innocuous stimuli, it was 3.6 (±2).
Temporal Clustering (Mechanical Pain)
Fuzzy clustering decomposed the whole-brain BOLD signal into 8 clusters (Fig. 1). The percentage of stimulus-locked brain voxels detected by this clustering method was 39%, while the GLM only showed 23% of the brain voxels as active, meaning that 58% of the voxels detected by the clustering method were characterized by a nontypical BOLD response. This result is in line with the findings recently reported by Gonzalez-Castillo et al. (2012), who described 68% of active voxels as having a nontypical BOLD response. These voxels are only partially detected by traditional fMRI analysis employing canonical HRF.
Of the 8 clusters, 2 were discarded because clearly composed of noise (Supplementary Fig. 2). The remaining 6 were considered as reflecting physiological signals.
The first cluster was composed of the primary somatosensory cortex, S1, the superior/middle temporal gyrus, the posterior insula, the mid-cingulate cortex, and the parietal and precuneal cortices. This cluster showed an early positive response to the stimulation closely resembling the 2-gamma BOLD responses (Fig. 2 and Supplementary Fig. 3).
The second cluster was composed of the anterior insula, the putamen, and the dorsal anterior cingulate cortex. Its activity was characterized by a wider and later positive response when compared with that of the first cluster.
The third cluster showed a fronto-parietal pattern of activity, including the dorso-lateral prefrontal cortex, frontal eye fields, and the premotor and temporo-polar cortices. Its temporal profile was characterized by a wide and sustained positive response.
The fourth cluster included portions of the default mode network (DMN; Raichle and Snyder 2007), and as expected, showed an early negative response (Hagmann et al. 2008). The fifth cluster principally involved primary, secondary, and associative visual areas. It showed a temporally widespread negative response (Supplementary Fig. 4).
The sixth cluster included voxels in the primary somatosensory, somatomotor, as well as premotor and somatosensory association cortices and was labeled “sensorimotor.” However, as the temporal profile of the sensorimotor cluster was characterized by a high similarity (Supplementary Fig. 5) with the temporal profile of the cluster which included the posterior insula, the sensorimotor cluster was considered an over-decomposition of this cluster.
The activations in the posterior insular cluster occurred earlier than those in the cingulo-insular and fronto-parietal clusters (Fig. 2) and, as previously anticipated, this cluster showed the canonical HRF. In contrast, the latter 2 clusters showed a more prolonged response, compared with the first one. The middle right panel shows a comparison between the temporal behavior of the GLM areas and the first cluster of painful stimulation. It can be observed that, although the profile of these 2 responses is similar, the GLM response is characterized by a more prolonged activation.
Bold Latency Mapping
Since we were principally interested in activations, only the 3 positive clusters (comprising the posterior insula, cluster 1; the cingulate and insular cortices, cluster 2; and the frontal and parietal cortices, cluster 3) were included in this analysis. An analytic investigation of the temporal characteristics of the 3 clusters confirmed the results of the event-related averages. The three clusters had a significantly different time-to-peak pattern, with the first cluster peaking before the cingulo-insular and fronto-parietal clusters. The onset latency of the posterior insular and cingulo-insular clusters was instead comparable, with only the onset latency of the “fronto-parietal” cluster showing a greater delay (Fig. 3, upper panel). The cingulo-insular and fronto-parietal clusters were found to have a significantly higher BOLD response (PSC). Although weak, the BOLD responses were considerably reliable. Indeed, the weakness of the signal was mainly due to the massive averaging as is clearly shown by the comparison with the single subject averaging (Supplementary Fig. 6).
We performed event-related averages and BLM, 2 methods leading to similar results, because we were interested in double-checking the temporal estimation of the hemodynamic response. Indeed, the techniques employed for parameter estimation in time-resolved fMRI suffer from the extreme spatial and temporal variability of the BOLD response. The comparison with a simple estimation as close as possible to the actual data, such as the event-related average, can provide a simple and effective means of double-checking the sophisticated BLM estimation.
Voxel-wise BLM (Fig. 3, lower panel) showed that the time-to-peak wavefront moved from the posterior insula, superior parietal, inferior frontal, and midcingulate cortices to the dorso-lateral prefrontal, anterior dorsal insular, supplementary motor, and anterior cingulate cortices.
The duration of the response was found to be greater for those areas characterized by a greater time to peak. The onset latency was more prolonged for the dorsal premotor, inferior parietal, and mesial premotor areas. The amplitude of the BOLD signal (PSC) was lower for the posterior insular component and higher for the cingulo-insular and fronto-parietal ones (for single-subject data and PSC analyses see also Supplementary Figs 7 and 8).
Temporal Similarity Analysis
The hierarchical clustering of the dissimilarity matrix allowed the creation of a hierarchical dendrogram of clusters. It can be observed that those areas characterized by a negative response are more similar to each other than those characterized by a positive response. In addition, the GLM activations are similar to those of the posterior insular cluster, and also to those of the cingulo-insular and fronto-parietal clusters (Fig. 4, upper panel), thus confirming that the GLM response captures all of these components.
The ISC of the 6 components and of the 2 artefactual components are shown in the left panel of Supplementary Figure 5. The values represent the percentage of the presence of the pattern for each subject (e.g. an ISC of 100 for one component means that such a pattern is present in 100% of the subjects). The ISC can therefore be regarded as an index of reliability. As can be seen, negative responses (deactivations) were characterized by lower reliability. The lower right panel shows a comparison between the posterior insular and the sensorimotor cluster response. The high similarity of the latter part of the response of these 2 clusters suggests that the 2 components may be part of a same over-decomposed component. Alternatively, the sensorimotor cluster may represent the deactivation of an effector network that is normally activated in case of threatening events, such as pain, for the purpose of avoiding harm. Given the intrinsically hierarchical nature of the brain clusters (Doucet et al. 2011), by imposing a different decomposition numerosity it may even be possible to find these 2 clusters to be distinct. Interestingly, the maximum negative PCS of the sensorimotor (sub)network is reached after the involvement of the other networks with putative sensory, saliency detection, and attentional purposes, supporting the interpretation of the sensorimotor network as possibly devoted to motor responses to the threatening stimulus.
BOLD Response Remodeling
To explore whether standard predictor-based fMRI analysis is able to disentangle the full spatio-temporal complexity of the brain response to a painful mechanical stimulus, we remodeled specifically tailored hemodynamic responses after the temporal clustering analysis. The remodeled response was derived from the mean response of each cluster and convolved with the stimulus presentation stick response to create specifically tailored predictors. The new predictors were employed for a new GLM random effect (RFX) analysis.
The results of the GLM employing the newly created regressors are shown in Figure 4, middle panel. Notably, while with the use of the new regressors we could disentangle the first (posterior insula) and third (fronto-parietal) clusters, this method was not able to isolate the second (cingulo-insular) cluster that always appeared collapsed together with the posterior insular cluster. This was probably because, in order to be correctly disentangled using this method, the activity in these areas, needs to be more separated in time.
When examining the specificity of the brain response to painful stimuli with Neurosynth (Fig. 4, lower panel), we found that only a minority of areas of the pain neuromatrix have a high specificity for pain: The highest specificity was observed in the dorsal posterior insula. Other regions such as the anterior insula, the temporal lobe, the cerebellum, the amygdala, pregenual, and subgenual cingulate, and prefrontal cortices showed a low specificity for painful stimulation.
The results of the control experiment showed that tactile stimulation activated a series of brain areas mainly comprising the primary and secondary somatosensory cortices. To evaluate the quality of the activations generated by the nonpainful control experiment, a GLM analysis was performed and the statistical maps of this analysis compared with a meta-analysis of touch experiments performed with Neurosynth. We found the results of our experiment to be very similar to those obtained with the meta-analysis (Supplementary Fig. 8) and to reveal a pattern of primary and secondary somatosensory cortices, motor, premotor, anterior cingulate, insular, prefrontal, inferior temporal, and parietal activations.
After the GLM analysis, a series of event-related averages were created for each positive cluster detected with the temporal clustering strategy. This analysis showed that, in the first cluster, the tactile stimulus elicited a temporal response similar to that of cluster 1 of the painful stimulation (see Fig. 5 and compare with the upper panel of Fig. 2). Conversely, cluster 2 was only minimally activated by the innocuous stimulus, the PSC elicited here was much lower and the temporal profile broader. The third cluster showed a similar response to that of the fronto-parietal cluster of the painful stimulation but with a reduced PSC in the first part of the response. Interestingly, whereas the first cluster showed a variance (expressed as standard error) that was very similar between mechanical painful and tactile stimuli, the other 2 clusters expressed a higher variance when involved in a nonpainful stimulation.
In this study, we found that, after painful stimulus is applied, a widespread pattern of brain areas characterized either by canonical HRF or nonconventional BOLD responses show statistically significant stimulus-locked responses. The areas showing a nonconventional BOLD response, and thus not easily detectable using traditional fMRI analysis techniques, were recently reported (Gonzalez-Castillo et al. 2012) to account for 68% of all the stimulus-locked areas. All the stimulus-locked activations (and deactivations) related to the painful stimulation used in this study accounted for 39% of the total brain volume, whereas the areas detected through GLM analysis accounted for just 23%. These stimulus-locked brain areas can be grouped into 5 clusters according to their time course of activation. Each of these clusters is comprised of different networks and has its specific onset latency, time to peak, duration latency and PSC. Our findings suggest that the areas that respond with stimulus-locked activation to painful stimuli are likely to reflect the activity of different networks, each having a different temporal behavior and, possibly, subserving different cognitive functions. Only a portion of those areas are reliably detected by canonical GLM fMRI analysis.
While numerous studies have investigated the spatial profile of activation in relation to painful stimuli (e.g. which areas generate a stimulus-locked response or activate following painful stimuli), fewer studies have been conducted on the temporal profile of the BOLD response. None, and this represents the main novelty of the present study, have investigated the stimulus-locked responses under the hypothesis of a global involvement of the brain in stimulus-related processing (Gonzalez-Castillo et al. 2012). Here, we used mechanical painful stimuli and analyzed the temporal profile of the BOLD response in a true data-driven fashion. Pin-prick stimuli applied with a hand-held pin-prick probe coactivate Aβ-, Aδ-, and C-fibers and therefore are less specific for nociception than laser stimuli, which selectively activate Aδ- and C-fibers (Bromm and Treede 1984; Treede 1995). However, notably, the surface of the mechanical probe is hardly relevant to single-unit responses of low threshold mechanoreceptors such as Aβ-fibers (Garnsworthy et al. 1988). In addition, different aspects of the pricking pain are coded differentially by Aδ and polymodal C-nociceptors (Garnsworthy et al. 1988; Tamura et al. 2004): Probe size and force are better reflected in the activity of the Aδ (Greenspan and McGillis 1994; Garell et al. 1996). Taken together, these observations allowed us to consider the stimulation as mainly reflecting input transmitted in Aδ primary sensory afferents. Indeed, previous studies have compared cortical activity following laser and pin-prick stimulations and reported highly similar spatial profiles (Baumgartner et al. 2010). Whereas heat pain has been studied much more in fMRI experiments, cortical responses to sharp (phasic) and tonic pressure, 2 kinds of mechanical pain, have also been characterized (Ringler et al. 2003). Phasic mechanical painful stimuli have been shown to activate the primary and secondary somatosensory cortices, the cingulate cortex, and the insular cortex (Ringler et al. 2003). Pain from tonic pressure (mediated by C-fiber mechanically insensitive afferents; Andrew and Greenspan 1999; Schmidt et al. 2000) elicited instead weaker responses that, however, also lasted during resting periods. Mechanical pain stimulation has been less widely used in fMRI studies, despite the fact that it is one of the most commonly used in animal studies.
The majority of the studies investigating the temporal profile of BOLD responses have focused on the BOLD time course in a limited number of regions (S1 and S2, Chen et al. 2002), and using a prespecified canonical HRF. As it has been shown (Aguirre et al. 1998; Duann et al. 2002; Neumann et al. 2003; Saad et al. 2003; Kannurpatti et al. 2010) that the hemodynamic properties of the BOLD signal may vary across brain regions, it may be argued that one possible confound of traditional fMRI analyses is that these rely on the assumption that cortical and subcortical responses to painful stimuli are constant across brain regions. Indeed, traditional GLM statistics assume a spatially invariant parametric model of the HRF (Friston et al. 1994; Smolders et al. 2007), and this may result in different sensitivity of different regions of the brain. Our method tries to overcome this issue by clustering all the possible shapes of BOLD responses into different groups of areas sharing a similar temporal profile. This method has been shown to permit an improved detection of stimulus-locked brain areas (Gonzalez-Castillo et al. 2012) and to better account for the spatial and temporal variabilities of the BOLD response (Smolders et al. 2007; Moulton et al. 2012). Many previous studies (Chen et al. 2002; Moulton et al. 2005) used long heat stimuli (≥9 s) and analyzed the datasets using a set of predictors convolved with a canonical HRF (however see Upadhyay et al. 2010; Moulton et al. 2012). By applying stimuli of short duration and with the use of unsupervised clustering techniques to investigate such spatio-temporal aspects of the BOLD brain responses to mechanical painful stimuli, we were able to identify 5 statistically significant clusters of active voxels sharing a common temporal profile; all these clusters represent a consistent part of the brain volume and partially overlap with known resting-state functional connectivity networks. Three of these components showed an increase in the BOLD signal after the stimulus was applied: cluster 1, which comprised the posterior insula and somatosensory cortices, cluster 2, made of the cingulate cortex and the insula, and cluster 3 formed by the frontal and parietal regions. These clusters were reliably detected in approximately 80% of participants. We suggest these 3 clusters may subserve different cognitive functions with cluster 1 reflecting the first readout of the characteristics of the mechanical stimulus and clusters 2 and 3 representing more of an attentional reorienting toward the potentially harmful stimulus (Legrain et al. 2011). The response of the first component resembled closely the early phase response recently observed by Moulton et al. (2012). However, their early phase response also included a part of our cingulo-insular cluster, probably as a consequence of failure by the GLM to discriminate between the temporal profiles of these 2 clusters. The differences between our results and those reported by Moulton and colleagues may also be ascribed to the kind of stimulus that was used. Indeed, the temporal profile of the BOLD response using mechanical painful stimuli has been studied far less (Lui et al. 2008). The cingulo-insular cluster was composed of areas that have been proposed as forming a “saliency detection network” (Seeley et al. 2007) and also of areas of the multimodal network for the detection of salience proposed by Downar et al. (2000, 2001, 2002). The fronto-parietal cluster was formed of areas described as part of a fronto-parietal control network (Vincent et al. 2008). Although anatomical definitions of such a network diverge (Corbetta and Shulman 2002; Shulman et al. 2002; Seeley et al. 2007; Corbetta et al. 2008; Vincent et al. 2008), the majority of authors converge on considering it to be purported to the reorienting of attention and cognitive resources toward salient stimuli for the organization of coherent behavior. Vincent et al. (2008) proposed that this network is devoted to the integration of information emerging from the external world and regarding internal states of the individual. While the time to peak of the first 2 components differed by about 1 s, there was no significant difference in their onset latency. In contrast, the third fronto-parietal cluster showed a greater magnitude of the BOLD signal and a later and more prolonged activation. This finding suggests that the clusters may subserve different and sequential cognitive functions, more related to the elaboration of the characteristics of the stimulus or more related to an attentional cognitive response to it. An alternative explanation to the interpretation of the second and third clusters as reflecting attention reorienting regards the possibility that their activations reflect the recruitment of C-fibers and may capture the activity of another sensory component, such as the first one. Indeed, studies comparing brain responses with first and second pain have consistently shown that, in line with the slower conduction velocity of C-fibers, cortical activity related to the activation of such nociceptors peaks approximately 2 s later than that related to the activation of Aδ-fibers (Ploner et al. 2002; Forss et al. 2005; Qiu et al. 2006; Veldhuijzen et al. 2009). However, the possibility that first and second pain are represented in different cortical networks is debated. Forss et al. (2005), using MEG, showed that first and second pain share a common cortical network. These considerations thus appear not to support the functional interpretation of the 2 clusters as reflecting “first” and “second” pain.
To understand how much this temporal profile reflected pain processing and how much nonspecific somatosensory processing, we performed a control experiment in which innocuous brushing was used. The results of the control experiment showed that, after tactile stimulation, it is possible to identify 3 clusters. The first, including the insular and somatosensory cortices, appeared as a functional equivalent of the first cluster of pain stimulation. This finding is in line with the results of previous studies (Upadhyay et al. 2010) and suggests that these areas perform a first, possibly sensory, evaluation of the somatosensory stimulus. Interestingly, however, clusters 2 and 3 obtained in the control experiment showed a reduced response, thus confirming the interpretation of the functional meaning of clusters 2 and 3 in pain as reflecting attention orienteering toward salient and potentially dangerous stimuli. What is more, Lui et al. (2008), by comparing cortical responses elicited by painful and nonpainful mechanical stimuli, observed a large spatial overlapping of the activated networks. Nevertheless, the authors reported higher activity in the mid-anterior insular cortex, mid-anterior cingulate cortex, and in the dorso-lateral prefrontal cortex following noxious stimuli. Our results support this view by showing that networks activating with a longer temporal delay and including insular and cingular regions do indeed show a greater PSC following painful rather than innocuous stimuli. Crucially, our findings support the existence of an intensity coding for somatosensory stimuli (Garcia-Larrea et al. 2010; Pomares et al. 2012), but they also suggest that it is unspecific (Mouraux et al. 2011) and thus make substantial contributions to the new conception of pain activations as related to salience detection or nonmerely pain processing (Iannetti and Mouraux 2010; Legrain et al. 2011, 2012; Mouraux et al. 2011; Cauda et al. 2012).
In response to painful stimuli, we identified the activity of 2 components composed of clustered voxels showing a decrease in the BOLD signal in response to the stimulus and then considered as “negative” or anticorrelated. Both components, comprising visual areas, and areas of the DMN network, may reflect the reallocation of cognitive resources to external stimuli (Sridharan et al. 2008). It has been suggested that when salient sensory stimuli have to be elaborated, the balancing of resources between the fronto-parietal and DM networks is altered in favor of the former (Sridharan et al. 2008; Vincent et al. 2008). Previous studies investigating the temporal profile of the BOLD response reported the presence of a double peak of activation specifically following noxious heat stimulation that was not present after innocuous heat or tactile stimulation (Chen et al. 2002; Moulton et al. 2005). In particular, the second, late response has been considered as specific for pain (see also Becerra et al. 2001; Moulton et al. 2005, 2012; Upadhyay et al. 2010). These authors suggested that late activations in S1, the anterior cingulate cortex, and the insula were specifically related to nociception as they were not detected after the application of nonnocieptive stimuli (Becerra et al. 2001; Moulton et al. 2005, 2012). In further support of the interpretation of the late response as specific for pain, in a recent study, Moulton et al. (2012) showed that ratings of perceived intensity coding, but not pain ratings, correlate with the activity of contralateral primary and secondary cortices, thus suggesting that heat intensity is likely to be coded in S1 and S2. Importantly, only the late phase of the hemodynamic response, which the authors identified in S2, correlated with the ratings.
Our results converge on the finding that the activation wavefront involved a later activity in the cingulate regions and anterior insula (Pomares et al. 2012). However, an alternative explanation for such activations may be related to attentional/control components (cingulo-insular and fronto-parietal components) rather than specifically to nociception. Another possible and simpler account for this difference relates to the stimulus used. Indeed, studies in which laser stimuli were applied have shown monophasic responses (Bornhovd et al. 2002; Buchel et al. 2002).
As already revealed by the ROI-wise BLM analysis, the only component having a later onset is the fronto-parietal one. The fronto-parietal component may reflect the redirection of attention toward salient stimuli and reflect the cortical underpinnings for the integration of external and internal information. This view is supported by electroencephalography (EEG) studies investigating responses to nociceptive stimuli. These studies have shown that the first responses observed after the presentation of a nociceptive stimulus (such as the N1, Bromm and Treede 1991) reflect the elaboration of stimulus-related characteristics (Mouraux and Iannetti 2009; Gallace et al. 2011). In contrast, later responses to the nociceptive stimulus (such as the N2 and especially the P2) are considered to represent multimodal responses related to attentional reorientering (Legrain et al. 2002, 2003; Mouraux and Iannetti 2009; Valentini et al. 2011; Torta et al. 2012), or awareness of the stimulus (Lee et al. 2009; Gallace et al. 2011; Torta et al. 2012).
The results of the temporal analysis emphasized the differences between classical GLM and model-free temporal clustering of the BOLD signal. GLM deactivations appeared to be closer, in their temporal profile, to the DMN component than to the visual component. In addition, these analyses emphasized the similarity between GLM activations and the first posterior insular cluster, rather than between GLM activation and the cingulo-insular and fronto-parietal clusters. This is important as it shows that, by employing ad hoc-modeled BOLD responses, the traditional predictor-related fMRI analysis method can disentangle the posterior insular and the fronto-parietal clusters, although it was not able to isolate the cingulo-insular cluster that always appeared collapsed together with the posterior insular cluster. Indeed, temporal analysis showed the cingulo-insular component to be similar to the posterior insular one. This explains why the GLM fails in teasing apart these 2 components. Importantly, our data-driven method was able to reveal that 39% of brain voxels showed stimulus-locked, statistically significant activity, while GLM analysis only showed 23% of the voxels as active.
With regard to the validity of our results, since the vascular structure (Bandettini and Wong 1997) and vascular coupling (Huettel et al. 2004; Devonshire et al. 2012) are not constant over different brain regions, the latency and shape differences we detected may also be explained by a series of nonneuronal and noise-related biophysical causes. The reasons for such variability have yet to be fully understood. In this sense, the present results should be interpreted with particular caution. Nevertheless, our findings would be of great interest even if they reflected a mere hemodynamic variability of vascular or anatomical origins. First, this possibility would beg an important question about why these well-organized patterns (that very interestingly resemble some of the resting-state networks) can be detected using variable envelopes of vascular or anatomical origin. If the origin of these large-scale networks (that are detectable during a task (Smith et al. 2009; Laird et al. 2012) as well as during the resting state (Cauda et al. 2011)) is (at least in part) vascular, or anatomical, then our interpretation of the results needs a substantial revision. However, this possibility is unlikely, as a recent paper (Chang et al. 2008) showed that removing the hemodynamic delay from the resting-state networks leads to a slightly different but globally preserved pattern of connectivity. Moreover, if the spatial variability of the hemodynamic delay biased our data (leading to false positives or hemodynamic-driven clusters), the same problem should also have biased GLM-based fMRI analysis (in this case leading principally to false negatives, but also false positives, and nonexistent anticorrelations).
However, while we cannot prove that all the detected activations are driven by true neuronal activity (Logothetis 2008), these and the following series of considerations led us to believe that the clusters could not be explained solely on the basis of motion artifacts or other physiological confounds. Indeed, (1) Papers inspecting the brain distribution of the hemodynamic delay have described a spatial pattern that is very different to the networks detected in this study (Saad et al. 2001; Chen et al. 2002); (2) the clusters we observed are characterized by spatial patterns that differ from the vascular harboring structure; (3) since our method was able to reliably separate the regions that showed the canonical HRF responses and presented a pattern of activations of the BOLD response, it can be considered valid; (4) response remodeling confirmed the validity of the temporal clusterization results; (5) it is plausible that different vascular structures may lead to different onsets, amplitude, and time-to-peak of the BOLD responses, but it is much more difficult to explain why a different vasculature should lead to responses showing this spectrum of positive, negative, peaked, and sustained envelopes; (6) as discussed above, some EEG results support our data, (7) these data are also in line with the recent discovery that, during an active task, the majority of brain areas show a stimulus-locked brain activity (Gonzalez-Castillo et al. 2012); (8) finally and crucially, the results of the control experiment showed that the tactile stimulus specifically and differentially modified the response of the cingulo-insular and fronto-parietal clusters when compared with the painful stimulus. This is fully in line with the (much) lower salience of the tactile stimulation and confirms the behavioral interpretation of the networks detected in this paper. As far as the TR is concerned, it should be considered that although higher fields can lead to a more temporally resolved pulse sequence (Menon 2012), the differences in response and in latency detected here are in the order of seconds and can thus also be detected with a 2-s TR.
An interesting (although speculative) alternative explanation posits that the baseline cerebral blood flow and BOLD response can have a strong effect on the BOLD response elicited by a subsequent stimulus (Davis et al. 1998; Hoge et al. 1999; Kastrup et al. 1999; Kim et al. 1999; Li et al. 1999; Corfield et al. 2001). As suggested by several other studies (Sapir et al. 2005; Ploner et al. 2006, 2010; Boly et al. 2007; Seminowicz and Davis 2007; Hesselmann et al. 2010; Pichè et al. 2010; Sadaghiani et al. 2010; Lee et al. 2011), it can be hypothesized that in our paradigm an interaction between large-scale brain networks (such as, for example, the DMN and fronto-parietal networks) and BOLD responses elicited by the stimuli may have led, in some specific areas, to a stimulus-locked response incorporating some of the temporal characteristics of the large-scale networks. This has been confirmed by previous electrophysiology studies which found that ongoing fluctuating brain activity exerts a significant modulatory effect on subsequent stimulus-evoked brain activity (Arieli et al. 1996; Makeig et al. 2002). That is, we suggest the possibility of a phenomenon of temporal synchronization between continuously ongoing brain fluctuations (related to large-scale brain networks) and BOLD responses to external stimuli such a synchronization would give origin to an envelope constituted by the convolution between ongoing fluctuation and BOLD responses. In this view, in areas belonging to specific large-scale networks, the stimulus-related BOLD signal is (at least partially) convolved with the network-specific resting time course to generate a complex signal that incorporates some characteristics of both. Previous studies on sensory perception (Ebner and Armstrong-James 1990) demonstrated that, at the cortical level, the perception of a stimulus depends on mutual enforcement of exogenous stimuli and endogenous inputs. This is further confirmed by the evidence that, modulating the influence of endogenous activity over the stimulus-related activity using different anesthetics (halothane/α-chloralose) on anesthetized rats, the brain activations differ substantially (John and Prichep 2005). Temporal clustering techniques such as the one utilized here may thus capture the spatio-temporal pattern of these areas, where a partial convolution of large-scale network-related time courses and stimulus-locked envelopes generate nonconventional BOLD responses (see Supplementary Fig. 9 for a graphical representation of the phenomenon).
In conclusion, we have shown that extensive activity of multiple networks is engaged at different time latencies after a painful mechanical stimulus is presented. These findings require further confirmation, but aim to motivate research on a controversial topic, such as the temporal profile of BOLD responses, the variability of these response profiles, and the interaction between the stimulus-related BOLD response and the ongoing fluctuation of the large-scale brain networks.
Conflict of Interest: None declared.