Abstract
This study addresses one long-standing question of whether functional separations are preserved for somatosensory modalities of touch, heat, and cold nociception within primate primary somatosensory (S1) cortex. This information is critical for understanding how the nature of pain is represented in the primate brain. Using a combination of submillimeter-resolution fMRI and microelectrode local field potential (LFP) and spike recordings, we identified spatially segregated cortical zones for processing touch and nociceptive heat and cold stimuli in somatotopically appropriate areas 3a, 3b, 1, and 2 of S1 in male monkeys. The distances between zones were comparable (∼3.4 mm) across stimulus modalities (heat, cold, and tactile), indicating the existence of uniform, modality-specific modules. Stimulus-evoked LFP maps validated the fMRI maps in areas 3b and 1. Isolation of heat and cold nociceptive neurons from the fMRI zones confirmed the validity of using fMRI to probe nociceptive regions and circuits. Resting-state fMRI analysis revealed distinct intrinsic functional circuits among functionally related zones. We discovered distinct modular structures and networks for thermal nociception within S1 cortex, a finding that has significant implications for studying chronic pain syndromes and guiding the selection of neuromodulation targets for chronic pain management.
SIGNIFICANCE STATEMENT Primate S1 subregions contain discrete heat and cold nociceptive modules. Modules with the same properties exhibit strong functional connection. Nociceptive fMRI response coincides with LFP and spike activities of nociceptive neurons. Functional separation of heat and cold pain is retained within primate S1 cortex.
Introduction
In sensory systems, functional separation of different sensory modalities, such as color vision in visual systems or light touch in somatosensory systems, permits fast information processing while preserving functional specificity. Loss or disruption of functional specificity during information processing may contribute to various pathological pain symptoms, such as perceptually distinct heat or cold allodynia and burning phantom pain (Lorenz et al., 2002; Flor et al., 2006). As a stand-alone sensory modality, whether such a segregated organization is preserved for thermal pain in the primate brain is still debated. Distinct psychophysical features of the pain sensation evoked by noxious heat and cold stimuli support the existence of separated channels for the modality (or nature) of pain (Davis, 1998). We focused this study on the primary somatosensory (S1) cortex because this region is believed to encode the sensory features (e.g., nature, temporal, and spatial) of heat- and cold-elicited pain (Davis et al., 1998; Flor et al., 2006; Erpelding et al., 2012; Makin et al., 2013; Liu et al., 2015). Changes to neuronal activity or somatotopy have also been linked to phantom and other chronic pain conditions (Flor et al., 2006; Wrigley et al., 2009; Gustin et al., 2012; Kim et al., 2017).
The overall aim of this study is to determine whether the functional separations are preserved in the S1 cortex of primates for touch and thermal (cold and heat) pain. One effective way to address this question is to perform a combination of fMRI and electrophysiology studies under identical experimental conditions. Such an approach is advantageous and necessary in several respects. First, fMRI has revolutionized our understanding of pain since its inception by allowing the opportunity to investigate the underlying neural mechanisms of pain sensation. It permits simultaneous measurements across multiple macroscale brain regions while the system is engaged in functions related to the perception and modulation of pain in healthy and disease states (Apkarian et al., 2005; Apkarian, 2010). The fundamental tenet of these imaging studies is that fMRI blood oxygenation level-dependent (BOLD) signals change in parallel with neural activity (Logothetis et al., 2001). Recent evidence, however, indicates that the relationships between fMRI BOLD and electrophysiological signals indeed vary for different brain regions or task conditions, perhaps due to differences in the functional organization of neurons within a specific imaging volume and/or differences in the engagement of neurons during each task (Rees et al., 2000; Logothetis et al., 2001; Mukamel et al., 2005; Shih et al., 2009; Bartolo et al., 2011). In this context, a better understanding of the neuronal constituents underlying MRI signals with respect to pain processing is essential for the full appreciation of their clinical and behavioral implications needed for building an organization-based mathematical model for reliably detecting pain-related fMRI signal changes in humans (Eklund et al., 2016, 2017; Cox et al., 2017; Kessler et al., 2017) and for selecting the intervention targets for neuromodulation with higher precision (Avenanti et al., 2005; Lee et al., 2017). Such multifaceted and mesoscale level (i.e., submillimeter to several millimeters) information is often impossible to obtain in human subjects for both technical and ethical reasons, but it can be acquired from animal studies. Nonhuman primates are an ideal model because their brain closely resembles the human brain in both structure and function, particularly in primary sensory areas (Hutchison and Everling, 2012; Chen et al., 2017; Shi et al., 2017).
In this study, we used a combination of submillimeter-resolution BOLD fMRI at 9.4 T, a single microelectrode, and 98-channel Utah array electrophysiology. We specifically (1) examined the spatial relationships of fMRI responses to innocuous tactile and nociceptive heat and cold stimuli; (2) compared the spatial correspondence between fMRI activations and local field potential (LFP) maps acquired with a 98-channel Utah array; (3) validated fMRI signal changes with spiking activity of a single nociceptive neuron; and (4) delineated the local, intrinsic, functional circuits of innocuous tactile, nociceptive heat and cold regions using seed-based resting-state fMRI (rs-fMRI) signal analysis. Here we provide fMRI, LFP, and single-unit spiking activity evidence supporting the presence of separated modules and distinct mesoscale circuits for processing innocuous tactile and nociceptive heat and cold information within primate S1 cortex.
Materials and Methods
Animal preparation
Seven male adult squirrel monkeys (Saimiri sciureus; SM-BK, SM-BW, SM-C, SM-H, SM-O, SM-R, and SM-V) were included in this study. All subjects underwent multiple fMRI scans and microelectrode electrophysiology mapping sessions. Two 7 × 7 (98) channel microelectrode arrays were used to obtain data in three monkeys (SM-O, SM-V, and SM-R). For fMRI experiments, each animal was initially sedated with ketamine hydrochloride (10 mg/kg)/atropine (0.05 mg/kg), intubated, and then maintained with mechanical ventilation and isoflurane anesthesia (0.5–1.1%) delivered in a 30:70 O2/N2O mixture. During scans, each animal was placed in a custom-designed MR cradle with head secured by ear and head bars. Physiology was maintained in a stable condition with a constant anesthetic delivery of isoflurane (between 0.7% and 0.8%). Vital signs, including peripheral oxygen saturation and heart rate (Nonin), end-tidal CO2 (22–26 mmHg; SurgiVet), and respiratory pattern (Small Animal Instruments), were monitored and recorded. Rectal temperature was monitored (Small Animal Instruments) and maintained between 37.5°C and 38.5°C using a circulating water blanket (Gaymar Industries). Intravenous administration of 2.5% dextrose in saline (3 ml/kg/h) was given throughout the imaging session to prevent dehydration and provide caloric energy.
For electrophysiology mapping experiments, the head of each animal was stabilized in a stereotaxic frame and was maintained under the same physiological conditions as fMRI scans. A round piece of skull was removed to expose the central and lateral sulci for microelectrode mapping. Cortex was stabilized and protected with 4% agar or silicone oil. Detailed procedures can be found in our previous publications (Chen et al., 2009; Wilson et al., 2016). All procedures were performed under aseptic conditions and were approved by the Institutional Animal Care and Use Committee of Vanderbilt University.
Stimulus Protocol
Thermal stimuli.
Under anesthesia, fingers were secured by gluing small pegs to the fingernails and fixing these pegs firmly in Plasticine, leaving the glabrous surfaces available for stimulation. Distal finger pads of digit 2 (D2 and D3) were stimulated with an Advanced Thermal Solutions thermal probe (16 × 16 or 30 × 30 mm2; ramp rate, 8°C/s; Medoc). To map nociceptive cold- and heat-evoked fMRI responses, we alternated blocks of baseline temperature (32°C, 30 s in duration) with nine blocks of either nociceptive cold (4°C and 7°C) or heat (47.5°C, 21 s in durations) within each imaging run. To quantify the cold temperature-dependent fMRI signal changes, we randomly alternated blocks of three temperatures (4, 7, and 15°C) with baseline (32 °C) blocks. Each temperature block was repeated five times in one run. Typically, within one imaging session (day), multiple imaging runs (3–9) were collected for each imaging paradigm (e.g., single or mixed temperatures). During fMRI scans, the thermode remained in contact with the skin during temperature changes.
Tactile stimuli.
Innocuous vibrotactile stimuli were provided by a rounded plastic probe with a diameter of 2 mm connected to a piezoelectric device (Noliac). Piezos were driven by Grass S48 square wave stimulator (Natus Neurology, Natus Medical) to indent digits in 0.43 mm vertical displacement at a rate of 8 Hz (pulse duration, 20 ms). Vibrotactile stimuli were delivered in 30 s on/off blocks to individual or a combination of distal finger pads (D1–D5). During stimulus off blocks, the probe was in light touch with the skin. Typically, each stimulus condition (e.g., D2) was repeated seven times within a single fMRI imaging session. Different digits were stimulated in a random sequence and typically were repeated 3–9 times (runs) within one imaging session.
fMRI data acquisition
All MRI scans were performed on a 9.4 T, 21 cm narrow-bore Varian Inova Magnet (Agilent Technologies) using a 3-cm-diameter surface transmit–receive coil positioned over the S1 and S2 somatosensory cortices contralateral to the stimulated hand. To maximize the mapping power, we placed four 2-mm-thick oblique image slices in parallel to the brain surface where the central and lateral sulci are located (Fig. 1a). T2*-weighted gradient echo high-resolution (0.156 × 0.156 × 2 mm3) structural images [TR = 200 ms; echo time (TE) = 14 ms] were acquired to visualize the blood vessel pattern on the brain surface. T2*-weighted (BOLD) functional MRI images (in 0.55 × 0.55 × 2 mm3 resolution) were acquired using a gradient echoplanar imaging (EPI) sequence (TR = 1500 ms; TE = 19 ms).
fMRI data analysis
Preprocessing.
Functional MRI signals (stimulation and resting state) went through standard preprocessing steps of slice timing (3dTshift, AFNI; RRID:SCR_005927) and 3D motion correction (3dvolreg, AFNI), and were corrected for physiological noise (respiration and cardiac) using RETROICOR (Glover et al., 2000). The stimulus-evoked fMRI EPI data were temporally smoothed with a low-pass filter with cutoff frequency of 0.25 Hz (fslmaths, FSL; RRID:SCR_002823) and then spatially smoothed using an isotropic Gaussian filter kernel with a full-width at half-maximum (FWHM) of 0.8 mm (3dmerge, AFNI). For rs-fMRI EPI data, a bandpass filter at 0.009 and 0.08 Hz was applied (3dTproject, AFNI). The motion parameters were also used as regressors in a general linear model to reduce their contributions to the rs-fMRI signal. Functional EPI images were upsampled from a 64 × 64 to a 256 × 256 mm2 matrix and coregistered with corresponding T2*-weighted high-resolution anatomical images using a linear image registration tool (FLIRT, FSL, FMRIB) for display.
fMRI activation maps.
fMRI activation maps were created using a cross-correlation function between the signal time courses of each voxel and the boxcar predictor of the HRF convolved stimulus presentation paradigm (3dDeconvolve, AFNI). Activation was defined by voxels that exhibited significantly correlated BOLD signal changes [p < 0.01, false discovery rate (FDR) corrected] and were organized in a minimum of five upsampled continuous voxels (cluster size of 0.78 × 0.78 × 2 mm3). The same criteria were used for both thermal and tactile activation maps. Thresholded fMRI activation maps (with statistical t values, typically with a threshold of t = 2) were spatially interpolated and then superimposed on the corresponding high-resolution T2*-weighted anatomical images using a linear image registration tool (FLIRT, FSL, FMRIB) for display purposes.
Generation of activation frequency map.
To generalize an activation pattern and evaluate the across-run reliability of activation maps to cold nociceptive stimuli, we generated frequency maps by quantifying the frequency of activations detected at each voxel (for details, see Chen et al., 2011). Activation maps of a single run were first thresholded at a t value = 2 (p < 0.01, FDR corrected). Then, a binary value of 1 was assigned to activated voxels (p < 0.01) and 0 to nonactivated voxels (p > 0.01). The activation frequency was then calculated for each voxel. The frequency maps (thresholded at a 30% probability) were used for quantifying the sizes of activations and the spatial overlaps between nociceptive cold and heat activation clusters. The spatial separation of different activations was quantified by the cortical distance between the centers of two activation clusters.
Resting-state fMRI data analysis.
The selections of tactile, heat, and cold seeds were based on the stimulus-evoked activation maps with anatomically defined area 3b regions or adjacent regions. Each correlation coefficient map was generated using one voxel as a seed; a statistical threshold of r > 0.45 was selected for display (see Fig. 6b–d; 3dfim+, AFNI). To investigate the relationships among different regions within tactile, heat, and cold resting-state functional connectivity (rs-FC) networks, inter-regional correlations were then calculated. Four S1 subregions (areas 3a, 3b, 1, and 2) corresponding to each correlation coefficient map were chosen for interareal correlation coefficient analysis. Selection of the ROI location of the face region in area 3b was estimated by the known somatotopic map and electrophysiology map for each animal. Pearson correlations among the five main cortical regions within these three networks were calculated on a single-voxel basis. A 13 × 13 region–region correlation matrix was generated using the mean correlation coefficient value of each region–region pair. Averaged inter-regional correlation coefficients with their corresponding SEs were calculated across all of the nine resting-state sessions from four animals and are presented as matrix plots in Figure 6, i and j. The values of the correlation coefficients associated with each ROI pair were taken for future statistical analysis at the group level. The statistical significance of the differences between the correlation coefficient values of ROI pairs was determined using a one-way ANOVA followed by Tukey's post hoc test. A result of p < 0.05 was interpreted to be statistically significant (Fig. 6k).
Quantification of BOLD signal time courses
We extracted the BOLD signal time course from the peak voxel (with a maximal t value) in each of the areas 3a, 3b, 1, 2, and M1 to quantify the amplitudes and temporal profiles of the BOLD responses to stimuli. Measures obtained in each run were then averaged across runs and animals and were examined for statistical significance using a nonparametric Kruskal–Wallis one-way ANOVA test followed by Dunn's method for multiple comparisons. A p < 0.05 was considered statistically significant. ROI-based BOLD time course results are presented as the mean ± SEM unless stated otherwise.
Alignment of cross-session fMRI activation and electrophysiology (blood vessel) maps
To compare fMRI activation obtained in different imaging sessions, and to validate fMRI activation with digit representation maps defined by microelectrode electrophysiology, we have developed methods to coregister the following different types of images: structural MRI images across imaging sessions (MRI–MRI), digital blood vessel images to structural MRI image (electrophysiology–MRI), and histology slices to blood vessel images (histology–electrophysiology). For the first two types of coregistration, we identified corresponding anatomical and blood vessel landmarks in each structural image, such as the visible surface vessels (dark strips) and transcortical veins (dark spots; Fig. 1d, example). These coordinates were then put into a point-based registration algorithm (implemented in MATLAB (RRID:SCR_001622); for details, see Chen et al., 2007; Shi et al., 2017). The registration transformation between these two sets of coordinates was then applied to the fMRI activation image, thereby coregistering the fMRI activation map to the structural MRI image (Hill et al., 1991; Zhang et al., 2007, 2010; Lecoeur et al., 2011). The same procedures were applied to align electrophysiology–MRI maps. Both the surface and transcortical blood vessel features were easily identifiable and used for registration of MRI maps with surface blood vessel maps obtained later during microelectrode mapping/recording sessions. The transcortical veins shown on the T2*-weighted images (as black dots) as well as electrode lesions (when they were available) were aligned with the vessel marks on the tangential histology section to coregister MRI–electrophysiology–histology maps.
Localization of fMRI activations in different cortical areas
We used anatomical, electrophysiological, and histological information obtained in each animal to localize fMRI activation foci to areas 3a, 3b, 1, and 2 and M1 cortex. Identification of fMRI activation in S1 subareas and M1 were usually straightforward after alignments of electrophysiology and fMRI activation maps. The rich and unique pattern of surface and transcortical blood vessels in each animal provided landmarks for coregistering maps obtained with electrophysiology and MRI. We confirmed fMRI and electrophysiology mapping results with histological evidence of recording sites by referencing the tissue lesion marks produced by the electrode array (Fig. 1c). The last key validation of fMRI activation came from electrophysiological microelectrode mapping and recording studies (Fig. 2a,b, selected penetrations shown by colored dots) in each animal. Information about neuronal response properties, including receptive field sizes, preferred stimuli (e.g., brush or light touch), and somatotopic organization, were used to determine cortical representations of digits and face in areas 3a, 3b, 1, and 2 of S1 (Chen et al., 2001; Sur et al., 1982). The digits (D1–D5) are represented as a palm-to-palm pattern at the area 3b–1 border and a digit tip-to-tip pattern at the area 3a–3b and area 1–2 borders. The approximate distance between digit centers in areas 3b and 1 is ∼1.7 mm. We focused our mapping on all five digit regions in areas 3b and 1 in all subjects, as well as in areas 3a and 2 in some cases. For nociceptive stimuli-evoked fMRI activation foci detected at the interareal border regions, we did not purposely localize them into a single area (e.g., area 3b or area 2), but instead classified them as activation modules at the interareal border zones.
Electrophysiology
Generation of LFP activation maps.
Broadband electrical LFP signals were recorded with two 7 × 7 (98 channels in total) Utah arrays using a multichannel Cerebus Neural Signal Processor System (Black Rock Microsystems). LFP signals were sampled at 500 Hz and then bandpass filtered between 0.1 and 150 Hz by a bandpass cheby1 filter for quantification. A 60 Hz notch filter was also used to remove power frequency interference. Voltage changes were measured against the signals of one reference electrode within the array. For each digit stimulation, a total of six to eight trials was recorded. Within each trial, 10 stimulus trains (each 30 s in duration) were presented with 30 s resting periods between trains. The power spectrums of LFP data were calculated at two different time periods: resting state and stimulating state. Frequency spectrum (1–150 Hz) was computed by fast Fourier transform. The response magnitude was calculated as a rate of gamma band power (30–150 Hz) during stimulation versus rest periods. We computed the locations and spatial extents of LFP to different stimuli (nociceptive heat and cold, and innocuous touch) responses in areas 3b and 1 (see Fig. 4, examples).
Single-unit recording.
Before recording, manual palpation with light tapping, nociceptive heat and cold stimuli were used to map receptive fields representing the hand and digits. To better localize the nociceptive neurons, the heat and cold fMRI activation maps were used as a reference. The receptive fields of neurons isolated at each penetration site were measured and recorded. The LFP signals and multiunit spiking activity were recorded by 8 or 16 channel linear electrodes (V-Probe, Plexon) using a Plexon multichannel recording system. Single units were isolated using off-line spike sorting with principal component analysis (Offline Sorter, Plexon; RRID:SCR_000012). Peristimulus time histograms (PSTHs) were then generated in response to various stimuli for quantifying the latency of onset, the total duration of responses, and maximum amplitudes (NeuroExplorer, Nex Technologies; RRID:SCR_0011818). Spiking rate increases that are above one SD of resting-state firing rate are considered as being responsive to stimuli.
Results
Distinct fMRI activation patterns to innocuous tactile, nociceptive heat and cold stimuli
We first mapped fMRI activations to different modalities of thermal nociceptive stimuli (heat vs cold) and related nociceptive activations to touch responses in each individual animal. Nociceptive heat and cold stimuli elicited multiple fMRI activation foci across the S1 subareas of 3a, 3b, 1, and 2, as well as in M1 cortex, in all subjects (Fig. 1d,f,g,i, two examples), whereas innocuous tactile stimulation of the corresponding distal fingerpads of D2 and D3 evoked fewer and more focal activations in areas 3b and 1 (Fig. 1e,h). The zoom-in 3D plots of the cold, tactile, and heat activation zones in area 3b (Fig. 1j–l) illustrate their discrete patterns. M1 activation was detected predominantly in nociceptive conditions. The composite maps showed clear spatial separation of the innocuous tactile and nociceptive heat and cold activation zones (Fig. 2). Using the well established somatotopic representation of digits in this species as a reference, the majority of nociceptive cold and heat activation zones (Fig. 2c, blue and red outlines) were located predominantly at or near interareal borders (e.g., areas 3b–1 and areas 1–2) with little spatial overlap, and did not overlap with tactile zones (Fig. 2a,b, green outlines, c, green domains). As a validation, robust and abundant low-threshold tactile neurons with appropriate D3 receptive fields were isolated (Fig. 2a,b, green dots) from the zones where D3 tactile stimulation evoked strong fMRI activation (Fig. 2a,b, green fMRI activation outlines). Similarly, nociceptive heat and cold neurons were isolated from the thermal fMRI activation zones in area 3b (see Fig. 4 for details).
Spatial relationships of fMRI activations in M1 and S1 cortices in response to the innocuous tactile and nociceptive cold and heat stimulation in two representative monkeys (SM-R and SM-H). a, T2*-weighted anatomical coronal image shows the placement of one oblique imaging slice centered around central sulcus (CS) for fMRI data acquisition. b, Side view of the schematic squirrel monkey brain shows the location of the fMRI image field of view. c, Cytochrome oxides stain of tangentially cut tissue section. Insert: zoom-in view of the microelectrode marks (red sticks) on a more superficial tissue section. Dotted yellow lines indicate area 3a-3b and area 3b-1 borders. Yellow arrows indicate the hand-face border. d, f, g, i, Multirun activation frequency maps to nociceptive cold (4°C; d, g) and heat (47.5°C; f, i) stimulation of digits 2 and 3 (thresholded at >33% probability). Color bar, Frequency of activation; 3/5, three of five runs. Scale bar, 1 mm. e, h, Corresponding tactile activation maps in the two animals (thresholded at t = 3; p < 0.001, FDR corrected). Color bar, t value range. j–l, 3D views of the cold (j), touch (k), and heat (l) activation patterns in areas 3b and 1.
Spatial relationships of innocuous tactile and nociceptive cold and heat activations within S1 (areas 3a, 3b, 1, and 2) and M1 cortices in four representative monkeys (SM-R, SM-O, SM-BW, and SM-H). a, b, Overlays of nociceptive cold (4°C), heat (47.5°C, color patches), and tactile (green outlines) fMRI activations on the electrophysiological map in monkeys SM-R (top) and SM-H (bottom). Color dots represent the electrode penetration site and receptive field properties of neurons (see color code for each digit and face). c, Composite maps of cold (blue), heat (red), and tactile (green) fMRI activations in four monkeys.
Spatial profiles and relationships between nociceptive cold and heat activations within S1 cortex
We found that the sizes of activation zones were comparable to those of different types of stimuli, with little spatial overlap among them. The mean cortical territory within S1 cortex that responded to nociceptive cold and heat stimuli of D2 and D3 was ∼25.97 mm2 (∼5.1 × 5.1 mm2) in size, regardless of the specific area within S1. Nociceptive cold and heat regions occupied 54.41% and 45.59% of the total nociceptive area, respectively, with a little overlap (4.66%; Table 1). These quantifications support the qualitative observations shown in Figure 2. The mean ± SD interzone distances between activation zones of different modalities within area 3b were comparable (2.48 ± 0.87, 2.40 ± 1.21, and 3.30 ± 1.87 mm) for tactile and cold, tactile and heat, and cold and heat, respectively). The interzone distances between the activation zones of the same modality were relatively larger than the distances between different modalities within area 3b (Table 2).
Activation sizes to nociceptive cold vs heat stimuli in S1 cortex
Interarea distances
Cold temperature-dependent BOLD fMRI signal amplitude changes in S1 subareas
The time course of fMRI signals derived from cold fMRI activation zones showed a typical 3–4 s delay between signal increase and stimulus onset. The responses to 4°C and 7°C stimuli were sustained during the stimulus presentation periods (21 s) and lasted >10 s after the stimulus ramped back to a baseline of 32°C. The measures of time to peak varied across areas (15, 25.5, 27, 27, and 27 s for areas 3a, 3b, 1, 2, and M1, respectively). The graded BOLD signal increases to 7°C versus 4°C stimulation (Fig. 3a, blue and green curves) were evident in all areas. Signal changes evoked by innocuous 15°C cool stimuli were very weak in general and fluctuated around the baseline (red curves) in some areas (e.g., areas 3b and 2). Pooled nociceptive cold (4°C and 7°C) signal changes were significantly stronger (0.6–1.1%) than innocuous 15°C responses (−0.2% to 0.2%) across all areas (Fig. 3g). The trends of temperature-dependent signal decreases (4–15°C) are present in all areas except area 3a, even though the signal differences between the two levels of nociceptive cold (4°C and 7°C) were not statistically significant (Fig. 3b–f). The strongest responses to 4°C stimuli occurred in M1 and area 2 (Fig. 3h,k). The response magnitudes to 7°C were comparable across areas (Figs. 3i,l), whereas those to 15°C were marginal (Fig. 3j,m).
Cold temperature-dependent percentage BOLD signal changes in different cortical areas. a, Group mean time courses of BOLD signals to nociceptive cold (4°C and 7°C) and innocuous cool (15°C) stimuli in areas 3a, 3b, 1, 2, and M1. Color shades around the color lines represent the range of SEs. Light orange background blocks indicate the stimulus presentation period (21 s). b–f, Plots of the percentage BOLD signal as a function of cold stimulus intensity (°C) in areas 3a, 3b, 1, and 2 of S1 and M1. g, Comparison of response amplitudes (percentage signal changes) to noxious cold (4°C and 7°C) vs innocuous (15°C) stimuli across areas. *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001 (nonparametric Kruskal–Wallis test followed by Dunn's post-test). h–m, Comparisons of the mean BOLD time courses (h–j) and peak (mean ± SE) BOLD signal changes (k–m) to nociceptive cold (4°C and 7°C) and innocuous cool (15°C) stimuli across cortical areas.
Spatial agreement between stimulus-evoked fMRI and LFP responses in areas 3b and 1
To determine whether the distinct fMRI activation patterns in response to different stimuli are indicative of underlying neuronal population activities, we compared the 2D-LFP activations directly with the fMRI activations in response to identical innocuous tactile and nociceptive heat and cold stimuli in three monkeys. The electrode recordings (Fig. 4e, left and right array boxes) revealed an orderly pattern of LFP power increases to tactile stimulations of D2 and D3 in area 3b (Fig. 5b–f, D1–D5 activation map). In area 3b, the LFP response evoked by tactile stimulation of D2 was located medially to that of D3 (Fig. 4a,b). Cold stimulus-evoked LFP responses were located at a more medial and posterior location to heat responses and at the border regions of the recording field of view (Fig. 4e, light blue patches). In area 1, the LFP activations of different digits were less well organized, but their locations in general agree with fMRI activations. Mainly, overlapping D2 and D3 tactile responses were detected in the bottom of the LFP recording maps (Fig. 4a,b, area 1 array on the right side). The nociceptive heat and cold stimulus-evoked LFP activities were detected at the left portions of the area 1 array, and centered at slightly different locations (Fig. 4, compare c, d, the right two area 1 arrays). Overlays of the heat, cold, and tactile LFP response maps show good spatial correspondence between fMRI and LFP responses (Fig. 4e). In supporting the robustness of 2D LFP recordings, Figure 5b–f shows the orderly organized D1–D5 representations in both areas 3b and 1.
The 98-channel (two 7 × 7 arrays) Utah array mapping of LFP responses to innocuous tactile (8 Hz) and nociceptive heat (47.5°C) and cold (4°C) stimulation in areas 3b and 1. a–d, Maps of LFP power increases in response to tactile stimulation of D2 (a) and D3 (b), and heat (c), and cold (d) stimulation of D2 and D3. Hand inserts show the stimulated digits. Color scale, Power (mV2). e, Overlay of fMRI (color patches) and LFP (color outlines) maps with respect to the electrode array (black dots) and digit representations identified by the receptive field properties of neurons (color dots). Scale bar, 1 mm.
Somatotopically organized LFP maps to tactile stimulation of individual digits (D1–D5). a, Spatial relationship between recording sites of Utah arrays and digit representations identified by single-microelectrode mapping and recording. Colored dots indicate the receptive fields of different digits. b–f, Maps of gamma power changes in response to tactile stimulation of individual digits. Color bar, Power ratio between stimulation versus baseline period. Scale bar, 1 mm.
Nociceptive heat and cold neurons were isolated in area 3b regions showing nociceptive heat or cold stimulus-evoked fMRI responses
We performed single microelectrode recordings to determine whether the nociceptive modality modules identified by fMRI contain heat- or cold-specific nociceptive neurons. In general, we found a high level of spontaneous activity at these recording sites. This finding drastically differs from the observation of very low spontaneous activity at the core tactile digit region in area 3b (Fig. 6, compare i, b). For example, at one heat fMRI activation cluster located at the border between areas 3b and 1 (Fig. 6a, red outline), we isolated a total of 80 single neurons (from two recoding sites: R1 and R2) and 10 of them (12.5%) were heat-sensitive (47.5°C) nociceptive neurons. The receptive fields of these heat nociceptive neurons were fairly large, and the boundaries were hard to define. PSTHs of four representative neurons show their response properties and location distribution to either distal pad or palm heat stimulation (Fig. 6b,f). At the R1 recording site, neurons that responded weakly to tactile stimuli were also isolated (Fig. 6c, left two PSTH plots). The right two PSTH plots in Figure 6c show the spontaneous firing of units isolated at the R2 recording site, which exhibited activities have no phase relationship to the stimulus presentation. Across animals, we found that the heat fMRI activation zones contained either heat only or a mix of heat and tactile neurons (Fig. 6i, schematic illustration). The table in Figure 6h summarizes the proportion of different modalities of neurons identified at heat versus cold fMRI activation zones. Similarly, at one zone near the border between areas 3a and 3b showing nociceptive cold-evoked fMRI activation (Fig. 6d, blue outlines and yellow dot), eight nociceptive neurons that responded to cold stimuli presented on either the distal pad of D2 or the palm were isolated (Fig. 6e, two units). No heat-sensitive or tactile neurons were identified at this location. The bottom two green units in Figure 6e show two examples of isolated spontaneously firing units during heat stimulation.
Single-unit spiking activities to innocuous tactile and nociceptive heat and cold stimuli at fMRI activation clusters in area 3b in three representative monkeys (SM-H, SM-B, and SM-K). a, d, Composite fMRI activation maps of nociceptive heat (red outlines), nociceptive cold (blue outlines), and innocuous tactile (green outlines) stimuli. Yellow dots labeled with R1 and R2 indicate the microelectrode recording sites. b, c, PSTH of isolated nociceptive heat (b) or tactile (c) units from two penetrations within the heat fMRI activation cluster (red outline) in area 3b of SM-H. The stimulation sites are illustrated on the insert hands for each unit. The shapes of isolated signal unit spiking are shown as inserts. Green horizontal lines indicate the mean firing rate (spikes/bin) at 95% confidence level. e, Nociceptive heat and cold units isolated from the cold fMRI activation cluster in area 3b of SM-B. f, Control microelectrode recordings (color dots) from the tactile stimulus-evoked fMRI activation cluster. g, PSTH of representative low-threshold tactile units in response to probe indentations. h, A table summarizing the total number of isolated heat, cold, and tactile stimuli-sensitive and non-stimulus-locked spontaneous unit activities at different fMRI activation clusters (heat vs cold). i, A schematic illustration of the relationships of the nociceptive heat and cold units isolated (color dots), the stimulation locations (thenar or digits), and their corresponding fMRI activation clusters (heat or cold).
The temporal characteristics of nociceptive neurons also differed significantly from those of classic tactile neurons in tactile zones of area 3b. The firing latencies of nociceptive neurons were generally varied and delayed, ranging from 0.5 to 10 s, with durations lasting from 5 s to beyond the 21 s duration of the stimuli (Fig. 6b). Tactile neurons isolated at the core digit region in area 3b (Fig. 6f, color dots); however, exhibited strong and transit firing activities at the stimulus onset and offset (Fig. 6g, examples). The temporal firing properties were not different between nociceptive heat and cold neurons, but were significantly different between tactile and cold nociceptive neurons in all four measures of time to peak, amplitude, duration of responses, and ratio of firing/stimulus duration; tactile and heat nociceptive neurons were different across three measures (Table 3). The firing rates of heat and cold nociceptive neurons increased and peaked slowly after stimulus onset. The relative ratios of firing/stimulus duration were close to 1.0 for heat (0.83) and cold (1.04) neurons, whereas that of tactile neurons was ∼0.5.
Temporal characteristics of spiking activity
Distinct mesoscale intrinsic functional circuits for nociceptive heat and cold within S1 cortex
Building upon our previous observations that highly associated cortical zones (e.g., D1 touch zones in areas 3b and 1) exhibit strong rs-FC (Wang et al., 2013; Wilson et al., 2016), here we examined whether functionally distinct and spatially separated innocuous tactile, nociceptive heat and cold modules form segregated and functionally specific mesoscale circuits. To illustrate the distinct rs-FC patterns of each modality zone, we placed three seeds at the voxel showing the strongest responses to nociceptive heat, cold, or tactile stimulation in area 3b (Fig. 7a, color dots) and performed a voxelwise functional correlation (i.e., FC) analysis. The overall correlation patterns of each seed revealed two features. First, the overall functional connectivity patterns for heat, cold, and touch seeds differed markedly (Fig. 7e, color patches) with limited overlap between them. Second, the rs-FC maps of areas 3b heat and cold modules overlapped locally (at the area3b/1 seed region) and to a varying degree in other areas (e.g., areas 3a or 2) to their corresponding fMRI stimulus-evoked activation zones (Fig. 7f–h, compare outlines and patches with the same color). To quantify to what degree the rs-FC strengths of functionally matched (e.g., heat to heat) versus functionally nonmatched (e.g., heat to cold) differ, we computed pairwise correlations across identified cortical zones by using highly correlated voxels in each area for each modality (e.g., heat modules in areas 3a or 1, 2 for each modality, 12 total) and a few voxels at face regions as non-hand control seeds (Fig. 7e, color dots). 2D matrix plots of correlation strengths among all possible combinations of seed pairs revealed several high-correlation clusters within each functional modality (Fig. 7i, each boxed cluster with dotted black outlines). Pairwise correlation quantification at the group level found that rs-FC strengths between the modules with the same functionality (e.g., area 3b tactile to area 1 tactile) were significantly stronger (p < 0.05, one-way ANOVA followed by Tukey's post-test; Fig. 7k, three groups of color-shadowed boxes) than those with different functionality (Fig. 7k, first-column groups). As a control, the rs-FC of the pairs within the same functional modality (e.g., touch) in hand regions were also significantly stronger (p < 0.001) than those of somatotopically nonmatched hand–face pairs (Fig. 7i, last blue column and row, k, gray columns). Building upon the group pairwise correlation analysis results, we summarize the main findings of the present study with a schematic illustration in Figure 8.
Distinct seed-based resting-state functional connectivity patterns of heat, cold, and tactile fMRI activation foci in area 3b of SM-H. a, Three seeds (colored dots) are selected in the touch (green), heat (red), and cold (blue) fMRI activation (color outlines) clusters. b–d, Corresponding voxelwise correlation (functional connectivity) patterns for each seed (labeled with *; thresholded at r > 0.45). Color bars indicate the r value (correlation coefficient) range. e, Overlaid heat, cold, and tactile seeds functional connectivity maps (color patches). f–h, Overlaid fMRI activation (colored outline) and functional connectivity maps (color patches) for tactile (f), heat (g), and cold (h) seeds, respectively. i, j, 2D matrix plots of pairwise correlation maps (i) and corresponding SE maps (j) across all modality and control (face region) seeds. T, Tactile seed in area 3b; T-3a, tactile seed in area 3a; T-a1, tactile seed in area 1; T-a2, tactile seed in area 2; H, heat seed in area 3b; H-3a, heat seed in area 3a; H-a1, heat seed in area 1; H-a2, heat seed in area 2; C, cold seed in area 3b; C-3a, cold seed in area 3a; C-a1, cold seed in area 1; C-a2, cold seed in area 2; F, face control seed. Color bar, r value range. k, Box plots of the r values (correlation coefficient) between different seed pairs. *p < 0.05, ****p < 0.001 (one-way ANOVA followed by Tukey post-test).
Schematic illustration of the functional modular organizations and mesoscale circuits within S1 cortex. a, Side view of functional areas of a new world monkey and the location of S1 (areas 3a, 3b, 1 and 2) and M1 cortex. b, Spatially segregated modular organization of innocuous tactile and nociceptive heat- and cold-processing regions and their functionally distinct mesoscale networks. Solid color lines indicate robust connections, and dotted color lines represent weak connections. 3a: area 3a; 3b: area 3b; 1: area 1; 2: area 2.
Discussion
Spatially discrete modules for heat and cold nociceptive information processing within S1 cortex
Emerging evidence from human and animal studies supports the notion that psychophysically distinct pain sensations evoked by noxious heat, cold, and mechanical stimuli (Chen et al., 1996; Davis, 1998; Morin and Bushnell, 1998; Green, 2004) are mediated through functionally specific peripheral and spinal neurons in ascending pathways (Basbaum and Woolf, 1999; Woolf and Ma, 2007; Lolignier et al., 2016). Whether or to what degree the functional segregations are preserved for pain generally and for pain modality specifically (evoked by cold, heat, and mechanical noxious stimuli) in the primate brain is still debated. Clinical studies support the presence of pain modality preferred brain regions. Lesions limited to the human thalamic principle somatosensory nucleus (Kim et al., 2007) or infarction of insula alter cold and heat pain sensation differently (Birklein et al., 2005). Human functional imaging studies have linked different whole-brain activation patterns to heat- versus cold-elicited pain sensations (Casey et al., 1996). Given the multifunctional nature of many of the identifiable cortical regions in the human brain (e.g., S1 cortex), different theoretical models have been proposed to account for the multidimensional features of pain perception (Apkarian et al., 2005; Moayedi and Davis, 2013). One of the remaining key questions is whether the nature of pain (e.g., burning pain evoked by noxious heat) is processed by shared or separated cortical regions and/or circuits.
By applying a combination of submillimeter-resolution fMRI, electrophysiology, and histology methods in individual monkeys, we discovered that innocuous tactile, nociceptive heat and cold inputs that originate from the same skin location on the body (i.e., digits here) are processed by spatially discrete neuronal modules within S1 cortex (Fig. 8, schematic illustration). We interpret the activation zones as functional modules based on two findings. First, the average cortical distances between fMRI activation foci are comparable (∼3.4 mm) for all sensory modalities (heat, cold, and tactile). This distance led to the estimated 1.7 mm radius of each activation module, which is comparable to known interareal distances between core digit touch zones in S1 (Chen et al., 2007). The modular cortical structure permits efficient information processing while retaining functional specificity. Second, the reasonable spatial correspondence between fMRI and LFP responses during the processing of different thermal nociceptive inputs confirms that the nociceptive-stimulus-evoked fMRI signal changes are of neuronal origin, from both a population (indicated by LFP) and single-neuron perspective. To our knowledge, this is the first evidence supporting the close relationship between fMRI and underlying neuronal population signal changes during the nociceptive processing at the mesoscale. Last, our most recent study of the same S1 tactile modules demonstrated that the interareal differences in rs-FC measures at the modular (or columnar) level covary with those of spontaneous low-frequency LFP activity, and the local spatial distribution of rs-FC signals is in close spatial agreement with that of LFP (Wilson et al., 2016; Shi et al., 2017). Together, the overall organization features we observed support the existence of millimeter-sized, functionally discrete modules and thermal nociceptive modality-specific processing circuits within the S1 cortex.
Distinct mesoscale functional circuits for processing innocuous tactile and nociceptive heat and cold information within S1
Numerous studies, including our own, have demonstrated cortical regions that are engaged in similar brain functions or are anatomically connected often exhibit strong rs-FC across multiple spatial scales (Fox and Raichle, 2007; Wang et al., 2013). Functionally (e.g., heat) and somatotopically matched modules (e.g., digit to digit) are strongly interconnected with each other. Strong interareal rs-FC related to strong neuronal functional connectivity measured by the coherence of local field potentials in early somatosensory areas (Wilson et al., 2016; Wu et al., 2017). One step further, here we discovered that heat, cold, and tactile modules form functionally distinct circuits. The presence of widely distributed, spatially discrete, somatotopically appropriate functional modules during stimulation with nociceptive heat and cold supports the hypothesis that the functional segregations for cold versus heat thermal nociception are retained at the first cortical station of S1 cortex. The thalamic inputs to these different modality zones within S1 subregions remain to be determined. It is possible that nociceptive specific inputs may arise from the thalamic nucleus VMpo discovered by Craig et al. (1994).
Engagements of areas 3b and 1 in the representation of the modality and temporal features of thermal nociceptive inputs
From an evolutionary point of view, complex behaviors exhibited by primates, including humans, are accomplished by streamlining different information into more specialized cortical areas for information extraction and integration. Specialized information-processing regions and pathways are the fundamental organizing principle of sensory systems, which permits preservation of functional specificity and improves information integration efficacy. One such classic example is the specialization of four highly functionally related but distinct subareas of 3a, 3b, 1, and 2 in S1 for haptic function (Kaas, 1993). Cutaneous tactile functions among these four subareas, such as frequency and spatial discriminations, are performed by area 3b (a homologous region of S1 in rodents) and area 1, whereas areas 3a and 2 are responsible for processing inputs from deeper receptors to form proprioception. The slow temporal features of those isolated heat and cold nociceptive neurons suggest that areas 3b and 1 are likely engaged in encoding both the slow and fast components of thermal pain sensations (Fig. 5). The nociceptive neurons do not appear to be involved in the encoding of the precise locations of the stimuli because thermal nociceptive neurons exhibited much larger receptive fields than their neighboring low-threshold tactile neurons, a finding that is in line with previous observations (Kenshalo et al., 1988). This suggests that it is unlikely that the precise localization of nociceptive heat- or cold-evoked sensation on digits is encoded by nociceptive neurons themselves, but rather by their interacting with nearby low-threshold tactile neurons that encode spatial location with much higher precision. It is possible that populations of nociceptive neurons may work together to encode the location of painful stimuli. Future studies are needed to determine this possibility. Another interesting finding is the high level of tonic spiking activity of the nociceptive neurons in the nociceptive modules. This feature is remarkably different from the low tonic firing activity exhibited by tactile neurons in core touch zones. The functional relevance of this high tonic firing activity remains to be explored. One note is that anesthesia may sharpen the spatial separations of different functional modules. We have shown with functional optical imaging of intrinsic signal that anesthesia reduces the receptive field sizes of touch neurons and the spatial extent of cortex responds to peripheral stimuli (Chen et al., 2009). We believe that the modular structures revealed under light anesthesia represent the core information process units in the S1 cortex of monkeys.
Correspondences of nociceptive stimuli-evoked fMRI, LFP, and single-unit responses: implication for probing nociceptive processing with fMRI
We found that both fMRI and LFP signals are more sensitive than spiking activity in detecting thermal nociceptive stimulus-evoked cortical responses. Although the unit isolation efficacy is much improved when an fMRI activation map is used as a guide, the total number of isolated nociceptive neurons is still relatively low in comparison with the abundance of tactile neurons in adjacent modules. We attribute these findings to several possible factors. First, the high tonic firing activity in the nociceptive domains often masks the stimulus-evoked firing activity, particularly when the firing onset lags by several seconds. Second, since synaptic activity and integration are the predominant contributors to fMRI and LFP signals, it is possible that nociceptive inputs are processed by neurons that are organized differently (e.g., the degree of functional homogeneity of neurons) or by neurons that belong to different interlaminar microscale circuits and/or mesoscale circuits. Spiking and LFP activities were sampled primarily from the middle layers (3–4), while the fMRI signals were sampled from the entire depth of cortex. Based on the known differences in the information flow across cortical layers within and between cortical areas, we speculate that heat and cold thermal nociceptive inputs, as well as tactile ones, may be processed by distinct interlaminar circuitries within and beyond S1 cortex. Future studies focusing on cortical layer-dependent recordings will provide key information for this speculation.
Conclusion
Here we report on fMRI and electrophysiology (spiking and LFP) evidence supporting the existence of discrete ∼1.7 mm modules and mesoscale functional circuits for the representation of somatosensory modalities of innocuous touch and nociceptive heat and cold within primate S1 cortex. Isolation of heat- and cold-sensitive nociceptive neurons from the fMRI activation zones in border regions between areas 3b and 1 further support the validity of using fMRI signals to probe thermal nociceptive stimulus-evoked brain responses and functional circuits. An improved understanding of the relationships among fMRI, LFP, and spiking activities in the context of pain processing have important implications for human fMRI studies. An improved understanding of the mesoscale functional organization and circuits for nociceptive information processing is critical for developing novel statistical models that can improve the specificity and sensitivity in detecting pain-related fMRI activity and circuits in the human brain. These findings are significant for both the pain research and functional neuroimaging communities, particularly for understanding the principles of nociceptive processing in cortex and the relationship between hemodynamic BOLD signals and the electrical activity of neurons at the modular level.
Footnotes
The present study is supported by National Institutes of Health Grant R01-NS-069909 to L.M.C. We thank Dr. Feng Wang and Fuxue Xin for assistance on fMRI data collection, Chaohui Tang for technical support on animal preparation, and George H. Wilson III for language editing of the manuscript.
The authors declare no competing financial interests.
- Correspondence should be addressed to Dr. Li Min Chen, Associate Professor, Institute of Imaging Science, Departments of Radiology and Radiological Sciences and Psychology, Vanderbilt University, AA 1105 MCN, 1161 21st Avenue South, Nashville, TN 37232. limin.chen{at}vanderbilt.edu