Abstract
The anterior insula has repeatedly been linked to the experience of aversive stimuli, such as pain. Previously, we showed that the anterior insula is involved in the integration of pain intensity and its prior expectation. However, it is unclear whether this integration occurs by a pain-specific expectation or a more general expectation of an aversive event. To dissociate these possibilities, we conducted an experiment using painful stimuli and aversive pictures with three levels of aversiveness on human male volunteers. Stimuli were preceded by a probabilistic, combined modality and intensity cue in a full factorial design. Subjective ratings of pain intensity and skin conductance responses were best explained by a combination of actual pain intensity and expected pain intensity. In addition, using fMRI, we investigated the neuronal implementation of the integration of prior expectation and pain intensity. Similar to subjective ratings and autonomic responses, the dorsal anterior insula represented pain intensity and expectations. The ventral anterior insula additionally represented the absolute difference of the two terms (i.e., the prediction error). The posterior insula only represented pain intensity. Importantly, the pattern observed in the anterior insula was only present if the cued modality was correct (i.e., expect pain); in case of an incorrect modality cue (i.e., expect aversive picture), the ventral anterior insula simply represented pain intensity. The stimulus expectation and prediction error specificity in the ventral anterior insula indicates the integration of expectation with painful stimuli in this area. Importantly, this pattern cannot be explained by aversiveness.
SIGNIFICANCE STATEMENT The anterior insula has been shown to integrate pain intensity and their expectation. However, it is unclear whether this integration is pain-specific or related more generally to an aversive event. To address this, we combined painful stimuli and aversive pictures with three levels of aversiveness. The ventral anterior insula represented pain intensity, expectation, and their absolute difference (i.e., the prediction error). Importantly, this pattern was only observed if the cued modality was correct. In case of an incorrect modality cue, this area simply represented as pain intensity. The stimulus expectation and prediction error specificity in the ventral anterior insula indicates the integration of expectation with painful stimuli in this area. Importantly, this pattern cannot be explained by aversiveness.
Introduction
Accurate and fast perception of pain is crucial for an organism to avoid tissue damage. In many cases, sensory cues precede pain and allow the organism to generate an expectation or prediction. The influence of such expectations on pain perception has been established in basic cued pain paradigms (Koyama et al., 2005; Keltner et al., 2006; Wiech et al., 2008; Atlas et al., 2010) and in the context of placebo and nocebo effects (Tracey, 2010; Häuser et al., 2012). It has previously been speculated that predictive coding, which integrates predictions and violation of predictions (i.e., prediction errors), could be a suitable mechanism mediating the effect of expectations on pain (Büchel et al., 2014; Geuter et al., 2017). In this framework, expectations (i.e., predictions) and the mismatch between these expectations and the sensory inputs (i.e., prediction errors) are assumed to be combined to minimize future prediction errors by updating the expectations at multiple levels of the neuronal hierarchy. So far, this framework has mainly been applied to visual and auditory processing (Rao and Ballard, 1999; Sterzer et al., 2008; Egner et al., 2010; Okada et al., 2018) and has only recently been extended to pain perception (Geuter et al., 2017). However, pain is not unimodal but comprises other aspects, such as aversiveness, unpleasantness, and salience, which are shared by negative stimuli in other sensory modalities. Therefore, the observed expectation and prediction error signals in the anterior insula and the intensity representation in the posterior insula (Geuter et al., 2017) might not be related to pain per se but could be related to aversiveness in general. This alternative hypothesis can be tested by adding a second aversive stimulus modality. Consequently, we adapted our cued pain paradigm (Geuter et al., 2017) by adding aversive pictures as a second sensory modality. In both modalities, we presented aversive stimuli at three levels of intensity. Each stimulus was preceded by a combined modality and intensity cue. Cues and stimuli were combined in a full factorial design with factors of cue modality, cue intensity, stimulus modality, and stimulus intensity. This design allowed us to investigate whether expecting an aversive stimulus with a specific intensity only affects the neuronal activity in response to a stimulus in the cued modality or whether this influence also generalizes to a stimulus in the other modality.
Within the framework of predictive coding, the prediction error represents the “surprise” induced by a mismatch between the sensory signals encountered and those predicted. If the previously observed prediction error signal in the anterior insula is related to aversiveness in general, we should observe such a signal if the cued intensity differs from the perceived intensity, regardless of the cued modality: a prediction error for a high-intensity painful stimulus, when a low-intensity aversive picture or a low-intensity painful stimulus was expected. On the contrary, if the anterior insula does not represent aversiveness in general, we should only observe prediction errors if the intensity cue is incongruent with the experienced intensity, but only if the modality cue matches: if a high intensity painful stimulus is perceived, when a low-intensity pain stimulus is cued, but not when a low-intensity aversive picture cue has been presented. Therefore, investigating the modality cue specificity of expectations and prediction errors for pain can potentially shed light on the characteristic of the pain processing system and a possible implementation of predictive coding therein.
Materials and Methods
Participants
Thirty-two healthy human male volunteers (average age: 26 years; range: 19–38 years) successfully participated in the study and were paid as compensation for their participation. Exclusion criteria were neurological, psychiatric, dermatological diseases, pain conditions, current medication, or substance abuse. All volunteers gave their informed consent. The study was approved by the Ethics board of the Hamburg Medical Association.
Stimuli and task
The stimuli were either thermal or visual stimuli. Thermal stimulation was performed using a 30 × 30 mm2 Peltier thermode (CHEPS Pathway, Medoc) at three different intensities: warm (42°C), medium painful (46°C), and highly painful (48°C) to the left radial forearm at different skin patches in each experimental block to avoid sensitization. The stimulation pattern was chosen based on the distribution of average human pain thresholds to have three distinguishable levels of thermal stimulation. The visual stimuli were pictures chosen from the International Affective Picture System (IAPS) (Lang et al., 2005) database at three levels of aversiveness. The classification of IAPS pictures was performed manually to divide the reported valance spectrum into six separate groups, where three of them were assigned to the training set and the other three were assigned to the fMRI set. The pictures in the fMRI set had three levels of aversiveness of which the low aversive category had aversiveness values of 2.02 ± 0.05 (mean ± SE), the medium aversive category had aversiveness values of 4.06 ± 0.02 (mean ± SE), and the high aversive category had aversiveness values of 5.23 ± 0.01 (mean ± SE). Before each picture or heat stimulus, a visual cue was presented. The color of the cue (triangle) indicated the modality of the stimulus (yellow for picture and red for heat). The size of the cue indicated its intensity (small, medium, and large triangle for low, medium, and high intensities, respectively).
At the beginning of each trial (Fig. 1), the cue was presented for 500 ms, indicating the modality and intensity of the following stimulus. The cue validity for the stimulus modality was 70%. Cue validity with respect to intensity was 60% (for details of all cue stimulus combinations and frequency, see Fig. 1F). After a jittered blank period of 300–700 ms, volunteers received the visual or heat stimulus with a duration of 2 s. After stimulus presentation, volunteers had to rate the intensity of the aversiveness of the picture or the heat stimulus within a time window of 2 s on a 4 point rating scale (1 labeled as “neutral” and 4 as “strong”). A fixation cross was presented throughout the trial. In addition, we added 4 catch trials to each block, in which volunteers were asked to report the preceding cue within 8 s. No stimulation was given in these trials. Volunteers were rewarded with 50 cents for a correct answer in each catch trial. Finally, a variable intertrial interval (2.5 ± 0.5 s) followed each trial, during which only a fixation cross was displayed.
Experimental design and predicted response patterns for INT, INT+EXP, INT+PE, and INT+EXP+PE models. A, In the INT model, the response to the stimulus is only driven by its intensity. B, In the INT+EXP model, the stimulus response is the weighted sum of the stimulus intensity (left) and the stimulus expectation (middle). C, In the INT+PE model, the stimulus response is the weighted sum of the stimulus intensity (left) and the stimulus prediction error (middle). D, For the INT+EXP+PE model, the stimulus response is the weighted sum of the stimulus intensity (left), the stimulus expectation (second from left), and the prediction error (third from left). Right, Examples shown here assume an equal contribution of the expectation and the prediction error signals to the stimulus response. Different contributions of the expectation or prediction error can lead to a different pattern for the stimulus response. E, Trial structure: At the beginning of each trial, a central cue indicated the probability of the intensity and the modality of the upcoming stimulus. Then after a blank screen with variable duration, volunteers received the stimulus (2 s), followed by the rating (2 s). The intertrial intervals were randomly jittered between 2 and 3 s. Volunteers were asked to fixate a central-fixation-dot (except for the rating duration). F, Table illustrating the cue-stimulus contingencies.
Procedure
The experiment consisted of two sessions: a behavioral training session and an fMRI session. Before the training session, volunteers were informed about the procedure and gave their written informed consent. They were informed about the experimental paradigm and the nature of the cues and the contingencies between the cues and the stimuli to avoid learning effects during the main fMRI session. The training session contained up to three blocks and was terminated once volunteers had learned the cue regularities. The termination criterion was defined as detecting >75% catch trials and a significant discrimination (repeated-measures ANOVA [rmANOVA]: p < 0.05) of the three levels of pain and aversive picture intensities. The fMRI session consisted of four blocks. Each block consisted of 130 trials (including four catch trials) and lasted ∼15 min. The trial order within each block was pseudo-randomized, and the order of blocks was randomized across volunteers. During each block, fMRI signal and skin conductance responses (SCRs) were measured. A high-resolution anatomical image for each volunteer was acquired at the end of the fMRI session. The whole fMRI experiment lasted for ∼2 h.
Data acquisition
fMRI data were recorded on a Trio 3 tesla system equipped with a 32-channel head coil (Siemens). Forty transversal slices (voxel size 2 × 2 × 2 mm, 1 mm interslice gap) were obtained within each volume using a T2*-sensitive single-shot EPI sequence (TR = 1.24 s, TE = 26 ms, flip angle: 60°, FOV: 220 × 220 mm) with parallel imaging (GRAPPA, in-plane acceleration factor 2) (Griswold et al., 2002) and simultaneous multislice acquisitions (“multiband,” slice acceleration factor 2) (Feinberg et al., 2010; Moeller et al., 2010; Xu et al., 2013) as described previously (Setsompop et al., 2012). The corresponding image reconstruction algorithm was provided by the University of Minnesota Center for Magnetic Resonance Research. Slices were tilted ∼30° relative to the AC-PC line to improve coverage in the brainstem. Additionally, T1-weighted structural images (1 × 1 × 1 mm resolution) were acquired using an MPRAGE sequence (TR = 2300 ms, TE = 2.98 ms, flip-angle = 9°). Stimulus presentation, response recording, and timing relative to MR data acquisition were done using Psychophysics Toolbox 3 (http://www.psychtoolbox.org). Skin conductance was recorded using an EDA100C MRI system (Biopac Systems), amplified using an analog amplifier (MP150, Biopac Systems), and sampled at 100 Hz using a CED 1401 analog-digital converter (Cambridge Electronic Design). SCR electrodes were connected to the palm (thenar and hypothenar eminences) of the left hand. During the study, lights in the MRI room were dimmed, and luminance was kept constant across volunteers.
Data analyses
SCRs.
The phasic skin conductance drive locked to the onset of each stimulus was estimated from the raw SCRs in each trial using an established deconvolution method as implemented in ledalab (Benedek and Kaernbach, 2010). The resulting phasic responses were then separately averaged for different conditions within a temporal window from 2 to 5 s after stimulus onset for the painful stimuli and within a window from 0 to 3 s for the aversive pictures, which was selected to cover the peak of the averaged SCR signal in each modality across volunteers.
fMRI data preprocessing.
Functional imaging data were analyzed using MATLAB (version 2014b, The MathWorks) and SPM12 (Wellcome Trust Centre for Neuroimaging, London). The first five volumes of each functional block were discarded. The remaining images were spatially realigned to correct for motion artifacts and were then spatially normalized using DARTEL, a high-dimensional warping algorithm implemented in SPM (Ashburner, 2007). The functional images were spatially smoothed using a Gaussian kernel with an FWHM of 6 mm. First-level models included separate regressors using a boxcar with a duration of 2 s for each of the 36 experimental conditions: 6 stimuli (3 levels of pain + 3 levels of aversive images) × 3 intensity cues × 2 modality cues. Additional regressors with a duration of 2 s modeled the rating period. Another regressor (boxcar, duration 9 s) modeled the catch trials. In addition, the model also included six motion parameters estimated during realignment. All regressors were convolved with a canonical hemodynamic response function as implemented in SPM. The first-level contrast estimates for the 18 pain conditions were the basis for our computational models. Model estimation used a linear mixed-effect model as implemented in MATLAB (version 2014b, The MathWorks). Models were estimated for all gray-matter voxels. Gray-matter voxels were determined using the segmentation routine in SMP12 at the single volunteer level, then averaged across volunteers, and finally thresholded at an 80% gray matter probability.
Neurological pain signature (NPS) analyses.
The NPS is a multivariate pattern of brain activity with the aim of distinguishing experimental pain from other conditions (Wager et al., 2013; Krishnan et al., 2016). The NPS response mirrors subjective pain and has been suggested as a surrogate marker for heat pain intensity (Wager et al., 2013). We computed NPS expression values for each of the experimental conditions based on the regressors described above separately for each participant (Wager et al., 2011, 2013; Krishnan et al., 2016). The resulting NPS values were then the basis for our model comparison.
ROI masks and analysis.
According to the finding reported by Geuter et al. (2017), insular cortex is the area among pain-responsive regions representing expectation and prediction error signals. Based on these findings, we mainly focused our analysis on the insular cortex and its anatomically defined subdivisions as defined by the Destrieux atlas (Fischl et al., 2004; Destrieux et al., 2010). Using insular subdivisions, we tested our four models on the average parameter estimates extracted from voxels within these ROIs (for the detailed model comparison routine, see below). Even though the insular cortex is the focus of this study because of its key role in the integration of pain expectation as the internal model, and sensory signals as the physical input (Ostrowsky et al., 2002; Nieuwenhuys, 2012; Geuter et al., 2017), the processing of nociceptive information is not limited to the insular cortex. We therefore conducted additional ROI-based analyses in several other pain-related regions, including the parietal operculum (SII), postcentral sulcus (SI), frontal medial cortex, anterior cingulate (ACC), posterior cingulate (PCC), amygdala, hippocampus, thalamus, cerebellum, and periaqueductal gray (PAG). Anatomical masks were the same as used by Geuter et al. (2017), except for the PCC, frontal medial cortex, hippocampus, and cerebellum, which were generated by thresholding the probability map of the corresponding areas in the Harvard-Oxford Atlas (Desikan et al., 2006) at 50% probability (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/Atlases). Similar to insular subdivisions, we tested our four models on the average parameter estimates from voxels within these masks.
Modeling behavioral, SCR, NPS, and voxelwise fMRI data.
Crossing our three levels of stimulus intensity, three levels of intensity cues, and two levels of modality cues resulted in 18 (3 × 3 × 2) different pain conditions, in which nine stimuli were preceded by a valid modality cue and nine stimuli were preceded by an invalid modality cue. For these valid and invalid cue conditions, subjective ratings, SCR, NPS, and voxelwise fMRI data were fitted separately to four different models. All models included a stimulus intensity term. The first model (intensity model [INT]) only contained the stimulus intensity term as the physical characteristic (intensity) of the stimulus (Fig. 1A). The second model (expectation model [INT+EXP]) additionally comprised an expectation term represented by the intensity cue (Figure 1B). The third model (prediction error model [INT+PE]), comprised two terms: the prediction error (i.e., the absolute difference between intensity cue and stimulus intensity) and stimulus intensity (Fig. 1C). This model was inspired by several studies reporting reduced neuronal activity in response to expected versus unexpected stimuli (Egner et al., 2010; Jiang et al., 2013). The final full model (INT+EXP+PE) contained both the expectation and the prediction error term, as well as the intensity term. Conceptually, this model assumes that the response to a painful stimulus is a weighted sum of its intensity value, the prediction error as well as the expectation (Fig. 1D). This model was inspired by studies showing the existence of both expectation and prediction error signals in different brain regions (Egner et al., 2010; Geuter et al., 2017). Prediction errors were coded in absolute differences, motivated by previous studies in the visual and auditory system (den Ouden et al., 2010; Todorovic et al., 2011; Todorovic and de Lange, 2012; Boll et al., 2013; Geuter et al., 2017). Furthermore, to account for the uneven “perceptual distance” between three levels of intensity, we added a second-order intensity term to all models. Model comparisons were performed on behavioral ratings, autonomous (SCR) data, as well as voxelwise BOLD signals and the NPS score. Model fit and parameter estimation were conducted using the FITLME function in MATLAB (version 2014b, The MathWorks).
Model comparison and statistical analysis.
To find the best of four alternative models, we used a nested-model comparison in a backward selection fashion. We started with the INT+EXP+PE model as the most comprehensive model and compared it with the INT+EXP and INT+PE models using likelihood-ratio tests. The likelihood-ratio test constructs a test statistic using the log-likelihood objective function evaluated under the unrestricted model parameter estimates (INT+EXP+PE) and the restricted model parameter estimates (INT+EXP or INT+PE). The test statistic follows a χ2 distribution with 1 df. Comparing the unrestricted INT+EXP+PE model with the two alternative-restricted INT+EXP and INT+PE models resulted in two p values. If both p values were below our threshold (e.g., p < 0.05), then INT+EXP+PE (unrestricted model) was selected as the best model. Otherwise, the model associated with the biggest p value was discarded, and the fit to the other restricted model was compared with the fit to the INT model as the simplest model using likelihood-ratio test. If the ensuing p value was below the threshold, then the intermediate model (in our case either INT+EXP or INT+PE) was the winning model. Otherwise, the INT model as the simplest model was determined as the best model. In summary, model selection proceeded according to the following algorithm:
1: Fit the data to the full model (INT+EXP+PE).
2: Compare the model with the two simpler models (INT+EXP and INT+ PE).
3: Find the largest associated p value from the previous tests.
≫ If the largest p value is above threshold, go to Step 4.
≫ Otherwise, determine the full unrestricted model as the best model and stop.
4: Compare the model with smallest associated p value with the INT model.
≫ If the p value is above threshold, determine the INT model as the best model and stop.
≫ Otherwise, determine the intermediate unrestricted model as the best model and stop.
The criteria used to select each model as the winning model are also illustrated in Figure 2. In addition, we used Bayesian model weights to verify the model comparison results of the likelihood-ratio test. The Bayesian model weight wM for the model M is given by the following:
where BIC is the Bayesian Information Criterions computed using FITLME function in MATLAB (version 2014b, The MathWorks). Each wM is between 0 and 1, and the model with the biggest wM is the best model. Furthermore, given the computed BIC, the Bayes factor comparing two alternative model was obtained as follows:
Therefore, the ratio of any two wM is equal to the corresponding Bayes factor comparing those models. For behavioral, SCR, and NPS data, the p value threshold was set at 0.01. At the fMRI voxel level, we used an ROI (insula) based correction for multiple comparisons using false discovery rate (FDR) correction at p = 0.05. FDR correction was performed using the FDR_BH function in MATLAB (version 2014b, The MathWorks). Moreover, we also tested for main effects and interactions across stimulus intensity, intensity cue, as well as modality cue using rmANOVA as implemented in MATLAB (version 2014b, The MathWorks). We also conducted paired t test, where all tests were two-sided unless otherwise stated. All statistical tests were conducted in MATLAB (version 2014b, The MathWorks).
Flowchart of the model selection procedure. A, The EXP+PE+INT model will be considered the winning model if the associated p values of the likelihood ratio tests are <0.05. B, The EXP+INT is considered the best model, if the p value for the likelihood ratio test of EXP+PE+INT against PE+INT is <0.05, but the p value for EXP+PE+INT versus EXP+INT is >0.05. C, The PE +INT model is considered the best model, if the p value for the likelihood ratio test of EXP+PE+INT versus EXP+INT is <0.05, and the p value for EXP+PE+INT versus EXP+INT is >0.05. D, If both p values of the likelihood ratio tests are >0.05, the INT model is considered the best model.
Results
Behavioral results
To validate our intensity manipulation for heat pain and to verify the discriminability between different levels of pain, we first tested for the main effect of stimulus intensity (Fig. 3A). Results regarding the aversive pictures are beyond the scope of this report but are depicted in Figure 3A for the sake of comparison. An rmANOVA, based on the behavioral data, revealed significant effects for the three levels of pain (F(2,62) = 736.4, p = 6 × 10−44, rmANOVA). Similarly, the main effects of intensity cue (F(2,62) = 30.7, p = 5 × 10−5, rmANOVA) and modality cue (F(1,31) = 10.85, p = 0.0025, rmANOVA) were significant. The two-way interaction between stimulus intensity and intensity cue (F(4,124) = 7.09, p = 3 × 10−5, rmANOVA) and the interaction between modality cue and intensity cue (F(2,62) = 11.45, p = 0.0001, rmANOVA) were also significant (Table 1). The rmANOVA results regarding the aversive pictures are summarized in Table 1.
Behavioral, autonomous responses and NPS to thermal pain and aversive pictures. A, Color bars represent the subjective ratings of warm (38°), medium-painful (46°), and highly painful (48°) stimuli (left) and subjective ratings for aversiveness of IAPS pictures (right). B, SCRs to the pain and aversive pictures. C, NPS responses to painful stimuli and aversive pictures. Error bars indicate SEM across volunteers.
Main effect of stimulus intensity, intensity cue, modality cue, as well as their corresponding two-way and three-way interactions: subjective ratings, SCR, and NPS in both pain and picture conditions
In addition, we fitted the behavioral ratings to our four alternative models (INT, INT+EXP, INT+PE, and INT+EXP+PE). Backward model selection revealed that the INT+EXP model best explained the data in the conditions with a valid modality cue. In contrast, the INT model best explained the data in the invalid modality cue conditions. In addition to the likelihood-ratio test, BIC and Bayesian model weights were computed for each model to validate the consistency of our model selection routine. Comparing the four models in the valid modality cue conditions based on their corresponding BIC, INT+EXP+PE model showed the biggest weight (BIC = 144.2 and BICw = 0.6), even though the Bayes Factor comparing INT+EXP+PE and INT+EXP model was BFPC>EXP = 1.52, which indicates their difference is not substantial; and as the expectation model (INT+EXP) has less degree of freedom, it can be determined as the best model describing the behavioral rating data in the valid modality cue conditions (Lewandowsky and Farrell, 2010). For the invalid modality cue conditions, again the INT+EXP+PE model had the biggest weight (BIC = 144.2 and BICw = 0.29), but there was no substantial difference between INT+EXP+PE model and the other three models (BFPC>EXP = 1.18, BFPC>PE = 1.21, BFPC>IC = 1.34). Therefore, the INT model as the simplest model was identified as the best model describing the behavioral ratings in the invalid modality cue conditions. In summary, both model selection strategies identified the INT model as the best model describing the behavioral ratings in the invalid modality cue conditions and the INT+EXP model as the best model for the data in the conditions with a valid modality cue.
SCR
Analogous to the behavioral data, the main effects of stimulus intensity (F(2,62) = 53.2, p = 3 × 10−14, rmANOVA), intensity cue (F(2,62) = 8.48, p = 0.006, rmANOVA), and modality cue (F(1,31) = 12.71, p = 0.001, rmANOVA) were significant. Similar to the ratings, the two-way interaction between stimulus intensity and intensity cue (F(4,124) = 12.33, p = 2 × 10−8, rmANOVA) as well as the interaction between the modality cue and the intensity cue (F(2,62) = 8.4, p = 0.0006, rmANOVA) were significant (Table 1; Fig. 3B). The best model describing the SCR data were identified using likelihood-ratio test in both valid and invalid modality cue conditions. As in the behavioral ratings, INT+EXP was the winning model in explaining the SCR data in the valid cue conditions. The INT model was the best model in the invalid modality cue conditions. This result was verified using Bayesian model weights for the model comparisons for the valid and invalid modality cue conditions. For valid modality cue conditions, the INT+EXP model had the biggest BICw of 0.78, whereas the BICw of INT+EXP+PE model was 0.21, the BICw of INT+PE model was 4 × 10−4, and the INT model had a BICw of 0.002. For the invalid modality cue conditions, INT+PE had the biggest BICw of 0.47, which was only 1.25 times bigger than the BICw of the INT model with the value of 0.37. INT+EXP+PE and INT+EXP models had BICws of 0.08 and 0.06, respectively. Therefore, the INT model can be considered as the best model in the invalid modality cue conditions.
fMRI data
In a first step, we analyzed fMRI data from all conditions using the NPS (Wager et al., 2013). For the NPS data, an rmANOVA showed a significant effect for the three levels of pain (F(2,62) = 122.65, p = 3 × 10−22, rmANOVA). Neither the main effect of intensity cue nor modality cue was significant. The two-way interaction between modality cue and intensity cue was significant (F(2,62) = 5.1, p = 0.008, rmANOVA); however, the other two-way interactions were not significant (Fig. 3C; Table 1). As for the behavioral data, we fitted the NPS-derived scores to all our four models. Analogous to the behavioral ratings, the INT+EXP model best explained the data for the valid modality cue conditions using hierarchical model comparison based on the likelihood-ratio test. In the invalid modality cue conditions, the INT model was best in explaining the NPS data based on the nested likelihood-ratio test. This was confirmed by Bayesian model weights in the valid and invalid modality cue conditions. The INT+EXP model had the biggest BICw of 0.88 (with the Bayes Factor of 14.83 compared with the INT+EXP+PE model as the second wining model) for valid modality cue conditions. The Bayesian model weights for INT+EXP+PE, INT+PE, and INT models were 0.06, 0.004, and 0.05, respectively. In the invalid modality cue conditions, the INT model had the biggest BICw of 0.89 (with the Bayes Factor of 16.15 compared with the INT+PE model as the second wining model). The other alternative models of INT+EXP+PE, INT+EXP, and INT+PE had the BICw of 0.003, 0.053, and 0.055, respectively. As control, we also tested the NPS on picture stimuli (Table 1; Fig. 3C).
We next investigated the response patterns of voxels in subregions of insular cortex according to a surface atlas (Fischl et al., 2004; Destrieux et al., 2010). Regions in this atlas are defined according to sulcal depth (Fischl et al., 2004; Destrieux et al., 2010). To obtain an intuitive visualization of the cortical representation for each model, we color-coded surface vertices according to the best explanatory model (four different colors) on the inflated cortical surface. As shown in Figure 4, the spatial distribution of the insular subregions according to this atlas closely resembled the distribution of the winning models. Activity patterns in the anterior division of the circular insula (ventral part of anterior insula [vAI]) were best explained by the INT+EXP+PE model, activity patterns in the superior division of the circular insula (dorsal part of anterior Insula) were best explained by the INT+EXP model, and the inferior division of the circular insula showed activation patterns compatible with the INT model. To verify this, we averaged the contrast estimates in each anatomical subregion of the insula and performed a formal model comparison on this average activation pattern. For the valid modality cue conditions, the INT+EXP+PE model best accounted for the fMRI pattern in the anterior division of circular insula, the INT+EXP model best described the response pattern in the superior division of circular insula, and the INT model was the best model explaining fMRI activation patterns in the inferior division of circular insula (Fig. 4, inserted). For all invalid modality cue conditions, the INT model was the best model across all insular subdivisions. We also compared the individual contribution of the expectation and the prediction error terms in voxels best explained by the INT+EXP+PE model within the insula. For this purpose, we mapped the corresponding parameter estimates within the insular cortex on the inflated cortical surfaces (Fig. 5A,B). This analysis revealed that the prediction error parameter estimate and the expectation parameter estimate are not significantly different in both hemispheres. We also observed that the expectation weight was not significantly different in INT+EXP representative cluster compared with INT+EXP+PE voxels in both hemispheres (Fig. 5C).
Voxelwise representation of the best model explaining responses to thermal pain within insular cortex. A, Contrast estimates were fitted to four alternative models of INT, INT+EXP, INT+PE, and INT+EXP+PE (valid modality cue). The INT+EXP+PE model (red) best explains the data in the anterior division of the circular insula (vAI). Color coding is based on likelihood-ratio test model comparison (FDR-corrected, q < 0.05). q values are the adjusted p values found using an optimized FDR approach in voxels within the insula. The INT+EXP model (yellow) best explains the fMRI data in the superior division of the circular insula (dAI). Within the inferior division of circular insula, the INT model (blue) is the winning model. B, When the modality cue was invalid, the contrast estimates were best explained by the INT model (blue) in all regions. Error bars indicate SEM across volunteers.
Voxelwise map of the respective weights of the individual components of the INT+EXP+PE and INT+EXP models within insular cortex for the valid modality cue conditions. A, Weights for the expectation term are plotted for all INT+EXP+PE representative voxels in insular cortex (i.e., Fig. 4, red voxels). Bar graph represents the average of the expectation weight within voxels best explained by the INT+EXP+PE model. Error bars indicate SEM across volunteers. B, Voxelwise map for the weight of the prediction error term as well as its average for all voxels best explained by the INT+EXP+PE model. C, Map of the expectation term and its average across voxels best explained by the INT+EXP model (i.e., Fig. 4, yellow voxels).
In an exploratory fashion, we also performed formal model comparisons with the threshold of 0.01 outside the insular ROI (Fig. 6) on the surface and in all voxels, including subcortical areas. Here we additionally observed groups of voxels representing the INT+EXP model as the best model in the ACC, but only in valid modality cue conditions (Fig. 6A). Among the subcortical areas, the PAG showed activation patterns that were best explained by the INT+EXP model (Fig. 6C). As in cortical regions, this was only true if the modality cue was valid. Otherwise, intensity coding was identified as the best model (Fig. 6B,C).
Color map representing the best model for each voxel at an uncorrected threshold of 0.01 projected on flattened cortex and inflated brain surfaces. A, In addition to the insular clusters depicted in Figure 4, an INT+EXP cluster is also evident within the ACC and visual cortex. There is also a small cluster of the INT+PE model within the cingulate cortex (green). B, For the invalid modality cue, the INT model is the best model across all pain-responsive areas. C, Color-coded model selection for subcortical areas. In parts of the PAG, the INT+EXP model best described the data. For the invalid modality cue, the INT model best explains PAG responses.
In an explorative fashion, we also examined the averaged parameter estimates in other pain-related anatomically defined ROIs for the valid modality cue conditions. Based on rmANOVA, we found a significant main effect of stimulus intensity in all ROIs, including ACC, PCC, parietal operculum, postcentral, frontal medial cortex, hippocampus, amygdala, thalamus, cerebellum, as well as the PAG (for detailed results, see Table 2). Moreover, among all ROIs, a significant main effect of cue intensity was observed in the PAG (F(2,62) = 3.88, p = 0.02, rmANOVA). The two-way interaction between modality cue and intensity cue was significant in the ACC (LH: F(2,62) = 4.62, p = 0.01; RH: F(2,62) = 4.63, p = 0.01, rmANOVA). With regards to model selection, we used condition-specific parameter estimates for each voxel and subjected those parameter estimates to our formal model comparison. As for the behavioral data, the model best explaining the data were identified using hierarchical likelihood-ratio tests. The INT+EXP model was determined as the best model (for detailed results, see Fig. 7; Table 3) in the ACC and PAG. In all other regions, the INT model was the winning model.
Main effect of stimulus intensity, intensity cue, modality cue, as well as their corresponding two-way and three-way interactions for all ROIs, including the SII, SI, frontal medial cortex, ACC, PCC, amygdala, hippocampus, thalamus, cerebellum, and PAGa
Parameter estimates of ROIs. Mean of the parameter estimates (± SE) at three pain intensities of high, medium, and low in different intensity cue conditions for valid modality cue trials are plotted for left (LH) and right (RH) hemispheres (except for the midline PAG) in each ROI. The fit of the parameter estimates to the best model among the four nested models is also depicted. In the ACC and PAG, the best model was determined as the INT+EXP model. Otherwise, the INT model was the best explanatory model. Bayesian model weights (BICw) for all different models are illustrated in each ROI. Error bars indicate SEM across volunteers.
Bayesian model weights for all ROIs, including the SII, SI, frontal medial cortex, ACC, PCC, amygdala, hippocampus, thalamus, cerebellum, and PAGa
Discussion
Our data show a clear discriminability of different levels of pain intensity based on behavioral ratings, autonomic measures, and a neural compound marker (NPS). Furthermore, we observed a strong effect of the intensity cue on pain perception, in agreement with previous findings (Koyama et al., 2005; Keltner et al., 2006; Atlas et al., 2010; Atlas and Wager, 2012; Wiech, 2016). Using computational modeling, activity patterns in the posterior insula were best explained by pure intensity coding, whereas activity in the dorsal anterior insula (dAI) was best explained by a model comprising the expected intensity. Finally, activity in the vAI was best explained by a model comprising stimulus intensity, cued intensity, and the prediction error (i.e., the difference between perceived and cued intensity) (Geuter et al., 2017). Finally, we observed a significant two-way interaction between modality cue and intensity cue, indicating a modality specific effect of expectation (Langner et al., 2011). Using formal modeling, we could reveal that the contribution of expectation in dAI and expectation plus prediction errors in vAI were only present for valid modality cues. In cases of an invalid modality cue, insular regions all showed simple intensity coding patterns.
Insular cortex
The insular cortex is crucial for pain processing (Ostrowsky et al., 2002; for review, see Nieuwenhuys (2012). It is ideally suited to represent pain intensity and related features because of its structural and functional organization (Barrett and Simmons, 2015; Nomi et al., 2017). While the traditional assumption (Cerliani et al., 2012; Cloutman et al., 2012; Jakab et al., 2012; Dennis et al., 2014) refers to an anterior–posterior dichotomy within insular cortex, more recent neuroimaging studies suggest a tripartite organization (Deen et al., 2011; Nomi et al., 2017).
In this tripartite organization, the posterior insula with its direct spinothalamocortical projections (originating in lamina I dorsal horn neurons) (Craig, 2002; Dum et al., 2009) as well as its intense connectivity to the somatosensory cortices (Wiech et al., 2014) makes it an ideal candidate for encoding the intensity of nociceptive thermal stimuli. Such sensory coding within the PI has previously been observed for both tonic and phasic pain (Segerdahl et al., 2015; Geuter et al., 2017) and even innocuous thermal stimuli (Davis et al., 1998). Our observation of a pure intensity coding pattern in the posterior insula confirms these findings.
The dAI is functionally connected to a cognitive control network (Dosenbach et al., 2007) and has stronger connections to posterior regions of the cingulate cortex (Taylor et al., 2009). Therefore, the dAI is presumably involved in expectation and attentional mechanisms (Deen et al., 2011). Moreover, the dAI as a part of cortical limbic areas with an agranular laminar structure is considered to encode expectation signals (Barrett and Simmons, 2015; Chanes and Barrett, 2016). This is in accordance with our findings showing a representation of stimulus intensity and expectation (INT+EXP model) in the dAI. Previous studies have also identified a role of this area in prestimulus effects and attentional modulation (Ploner et al., 2010; Wiech et al., 2010).
In contrast, the vAI exhibits strong connectivity to the posterior insula, ACC, thalamus, amygdala, and basal ganglia, which are also contributing to pain processing (DuPont et al., 2003; Baliki et al., 2009; Mutschler et al., 2009; Uddin et al., 2014). Such a strong connection of the vAI to these limbic structures makes it an ideal mediator between areas representing pain intensity and expectation to evaluate the expectation against the sensory input and estimating a prediction error (Downar et al., 2002; Seymour et al., 2004; Geuter et al., 2017). Our observed representation of stimulus intensity, expectation, and prediction error (INT+EXP+PE model) in the vAI is in agreement with these findings.
Expectation and prediction error in the valid modality
Predictions and prediction errors are key ingredients for perceptual theories, such as predictive coding (Knill and Pouget, 2004), and theories of brain function, such as the free energy model, including active inference (Knill and Pouget, 2004; Friston, 2005, 2008). In the predictive coding framework, the brain is relying on its internal model of the physical reality but tries to make it an adequate model by updating it based on the incoming data. Consequently, the integration of “prediction errors” into this model will change it over time to better reflect the physical reality. These concepts can account for extra-classical-receptive-field effects (Rao and Ballard, 1999), repetition suppression (Summerfield et al., 2008; Todorovic and de Lange, 2012), illusory contours formation (Dura-Bernal et al., 2011), and mismatch negativity (Ylinen et al., 2016). Such a mechanism was also proposed as a model to account for expectation effects in pain (Geuter et al., 2017) and placebo hypoalgesia (Büchel et al., 2014; Grahl et al., 2018). According to predictive coding in the case of pain, we hypothesize the presence of sensory specific representations of the expectation and prediction error terms at the hierarchical level where the intensity of pain is encoded similar to other modalities (Rao and Ballard, 1999; Friston, 2005). Crucially, this predicts that these representations are restricted to the case in which the modality cue is valid. As our second stimulus category comprised aversive pictures, this allows to test whether aversiveness, which is shared by pain and picture stimuli, can account for insula activation.
Importantly, our data show that the prediction error signal in the anterior insula is not generally related to aversiveness, which is shared by the pictures and the painful stimuli, but is limited to the pain expectation. This is indicated by the fact that we could only observe the “expectation” effect in the valid modality cue conditions. In case the modality cue was invalid (when volunteers expected an aversive picture but received pain), the intensity cue had no effect on the response patterns in the insula. In other words, an expectation that shares an aversiveness component (i.e., aversive pictures) is not encoded in these areas. This is further corroborated by the absence of a prediction error in the insula in this case. As a consequence, our data suggest a more specific role for the vAI, representing expectations and prediction errors related to pain. This indicates that this area is part of the hierarchy where pain intensity integration occurs and can be seen in analogy to other stimulus modalities, in which feature-selective predictive coding within the fusiform face area in response to faces, and within the parahippocampal place area in response to houses (den Ouden et al., 2010; Egner et al., 2010; Jiang et al., 2013) has been shown. Moreover, the observed functional dissociation between insular subregions resonates with previous accounts of insular subspecialization for different features of aversive stimuli (Baliki et al., 2009; Deen et al., 2011; Nomi et al., 2017).
Signed versus unsigned prediction errors
Prediction errors could either be represented as absolute deviations (i.e., absolute deviation of outcome from expectation) or signed deviations (i.e., more or less pain than expected). Unfortunately, due to the design of our experiment, distinguishing between unsigned and signed prediction errors cannot be dissociated in our model. Our data indicate the representation of unsigned prediction error in the anterior insula. Similar observations of unsigned prediction errors have been made in the anterior insula, amygdala, and visual and auditory cortices (den Ouden et al., 2009, 2010; Todorovic et al., 2011; Todorovic and de Lange, 2012; Boll et al., 2013; Geuter et al., 2017). Unlike the signed prediction error, which can be used to guide approach and avoidance behavior, this type of error (i.e., absolute prediction error) determines the captured attention, which is related to stimulus saliency rather than its actual value (Kahnt et al., 2014).
In conclusion, the vAI exhibits all features necessary for the implementation of a predictive coding type integration of pain. Importantly, the functional organization revealed by our model comparison is remarkably similar to the structurally defined tripartite organization of the insular cortex (Fischl et al., 2004; Destrieux et al., 2010). Our results corroborate the previously observed gradient from stimulus intensity coding to perceptual coding from posterior to anterior insula (Geuter et al., 2017). This might also reflect a gradient from nociception to pain, where somatosensory information from posterior somatosensory regions is more and more integrated with cognitive factors such as expectation, which finally allows the estimation of prediction errors in the anterior insula.
Footnotes
This work was supported by European Research Council Advanced Grant ERC-2010-AdG_20100407 and Deutsche Forschungsgemeinschaft Grant SFB 936 Project A06. We thank the University of Minnesota Center for Magnetic Resonance Research for providing the image reconstruction algorithm for the simultaneous multislice acquisitions; Tor D. Wager for sharing the NPS pattern; and Jürgen Finsterbusch, Katrin Bergholz, Kathrin Wendt, Waldemar Schwarz, and Matthias Kerkemeyer for valuable assistance with data collection. S.F. thanks Bjoern Horing for help on the NPS analysis, discussions, and comments on the manuscript draft.
The authors declare no competing financial interests.
- Correspondence should be addressed to Dr. Sepideh Fazeli, Department of Systems Neuroscience, University Medical Center Hamburg-Eppendorf, Martinistrasse 52, Hamburg, Germany 20246. sfazeli{at}uke.de
