Alterations of distributed neuronal network oscillations during acute pain in freely-moving mice.

The experience of pain involves the activation of multiple brain areas. Pain-specific activity patterns within and between these local networks remain, however, largely unknown. We measured neuronal network oscillations in different relevant regions of the mouse brain during acute pain, induced by subcutaneous injection of capsaicin into the left hind paw. Field potentials were recorded from primary somatosensory cortex, anterior cingulate cortex (ACC), posterior insula, ventral posterolateral thalamic nucleus, parietal cortex, central nucleus of the amygdala and olfactory bulb. Analysis included power spectra of local signals as well as interregional coherences and cross-frequency coupling (CFC). Capsaicin injection caused hypersensitivity to mechanical stimuli for at least one hour. At the same time, CFC between low (1-12 Hz) and fast frequencies (80-120 Hz) was increased in the ACC, as well as interregional coherence of low frequency oscillations (< 30 Hz) between several networks. However, these changes were not significant anymore after multiple comparison corrections. Using a variable selection method (elastic net) and a logistic regression classifier, however, the pain state was reliably predicted by combining parameters of power and coherence from various regions. Distinction between capsaicin and saline injection was also possible when data were restricted to frequencies <30 Hz, as used in clinical electroencephalography (EEG). Our findings indicate that changes of distributed brain oscillations may provide a functional signature of acute pain or pain-related alterations in activity.

Introduction

Pain is by definition a self-reported, subjective state. On the other hand, it is clearly associated with physiological processes both in the peripheral and central nervous system. Deriving a neurophysiological signature of pain from characteristic neuronal activity patterns would be very desirable to enable reliable diagnosis and distinction of different subclasses of pain (Davis et al., 2017), assessment of pain in patients with impaired communication (Benoit et al., 2017; Hartley et al., 2017), longitudinal monitoring of therapies and development of biofeedback methods (Jensen et al., 2013). In the brain, pain goes along with activation of multiple areas such as the anterior cingulate cortex (Lieberman and Eisenberger, 2015), the somatosensory cortex (Gross et al., 2007; Dowman et al., 2008; Zhang et al., 2012), the insular cortex (Frot and Mauguière, 2003; Brooks and Tracey, 2007; Isnard et al., 2011; Segerdahl et al., 2015) and the amygdala (Simons et al., 2014; Corder et al., 2019), summarized as the ‘pain matrix’ or ‘pain network’ (Melzack, 1999). Due to the broad distribution of participating networks it is difficult to isolate the contribution of single regions to pain processing (Mouraux et al., 2011; Davis et al., 2015; Wager et al., 2016). It rather seems appropriate to approach pain by a multiregional approach. As a first step in this direction we measured network activity in a wide range of brain regions in mice with and without acute pain, analyzing coordinated activity patterns within and between different local networks.

We focused on state-dependent network oscillations as a promising candidate for a neuronal signature of pain. Pain-dependent changes in electrical activity have already been described in several brain regions (Ploner et al., 2017). Such oscillations occur synchronously in different, functionally linked areas and are associated with specific functional, behavioral or vigilance states (Engel and Singer, 2001; Buzsáki and Draguhn, 2004; Uhlhaas and Singer, 2006; Fries, 2009; Buzsáki and Wang, 2012; Singer, 2018). It is therefore well possible that activation of a distributed functional system is reflected by a specific pattern of neuronal network oscillations throughout the brain. In addition, network oscillations form higher-order patterns which enrich the parameter space for potential functional signatures. For example, low-frequency rhythms such as theta oscillations can modulate high-frequency oscillations like gamma (Canolty et al., 2006). Such cross-frequency coupling (CFC) is linked to specific cognitive activities (Tort et al., 2010) and may well contribute to the signature of complex cognitive-emotional states such as pain. Indeed, CFC has been suggested as a marker of Parkinson’s disease (de Hemptinne et al., 2013) and of animal models of Alzheimer’s disease (Zhang et al., 2016).

Here, we used multi-site recordings of field potentials in the mouse brain to search for a signature of pain which was induced by injection of capsaicin into the hind paw. The resulting ‘tonic pain’ sensation is slightly prolonged compared to the shortest forms of phasic pain (lasting milliseconds to seconds), allowing to measure altered oscillation patterns during seconds to minutes. It is, however, clearly distinct from chronic pain. Mouse models of tonic pain are, hence, a good starting point to search for an electrographic signature of acute pain. Oscillations were measured in seven brain regions and covered six different frequency bands. Analysis included power within each frequency band, interregional coherence and cross-frequency coupling. Additionally, respiration rate and motor activity were recorded to exclude confounding factors. We report that capsaicin-induced pain is correlated with a specific combination of electrographic parameters and, hence, with a specific distributed pattern of network oscillations throughout the mouse brain.

Materials and methods

Local field potential (LFP) activity was obtained from depth electrodes in five pain-related brain regions, namely the primary somatosensory cortex (S1), anterior cingulate cortex (ACC), ventral posterolateral thalamic nucleus (VPL), posterior insula (Ins) and central nucleus of the amygdala (AMYG) and in addition in the olfactory bulb (OB) and parietal cortex (PAC). Due to channel restriction in our setups AMYG, OB, and PAC were only recorded from the left hemisphere (ipsilateral to injection side). Measurements were taken previous to and following local injection of saline (control) or capsaicin (tonic pain) into the left hind paw. Changes of LFP activity in six frequency bands (1−4 Hz, 4−12 Hz, 12−30 Hz, 30−80 Hz, 80−120 Hz, 120−160 Hz) were analyzed using power spectrum density (PSD), cross-frequency coupling (CFC) and coherence analysis.

Ethics statement

The experiments have been approved by the Governmental Supervisory Panel on Animal Experiments of Baden Württemberg, Karlsruhe (G-115/14). All efforts were made to minimize the suffering of the animals as well as the number of animals used for these experiments. There was no alternative to in vivo studies due to the systemic scope of this research.

Animal housing

C57BL/6 N mice were delivered from Charles River at an age of 10–12 weeks. The animals were housed in groups of four until surgery with ad libitum access to water and food. After surgery the mice had to be housed separately to avoid damage of electrode implantations by other mice. Animal cages were kept in a ventilated Scantainer (Scanbur BK) with inverted 12 h day/night cycle (lights on from 8 pm to 8 a.m.). After completion of experiments the mice were killed during brain perfusion by an overdose of pentobarbital (see below).

Animal preparation

Sixteen animals (6 female, 10 male) were anesthetized with isoflurane in medical oxygen (4 % for induction and 1–1.5 % for maintenance, flow rate 1 l per min). For analgesia, buprenorphine injections were applied after induction of anesthesia and repeated 8 h after surgery with a dosage of 0.1 mg per kilogram body weight.

The animals were placed in a stereotaxic apparatus (Kopf Instruments) while the animal’s body temperature was controlled by a heating pad (ATC-2000, World Precision Instrument). After exposing the animal’s skull, 0.6 to 0.7 mm holes were drilled according to the stereotaxic coordinates (see Table 1 and Fig. 1) without damaging the dura mater. For epidural recordings, stainless steel watch screws (1 mm diameter, 3 mm length) were positioned above the parietal cortex (PAC). Two additional surface electrodes were placed above the cerebellum serving as ground and reference electrodes. Depth electrodes, made from pairs of insulated tungsten wire (50 μm diameter; California Fine Wire, USA), were glued together and implanted in the remaining brain regions according to stereotaxic coordinates (Paxinos and Franklin, 2004, see Table 1) For monitoring respiration, thermal sensors (80 μm diameter, Omega Engineering Ltd, Germany) were implanted into the nasal cavity as described previously (Jessberger et al., 2016). All electrodes were connected to gold-covered pins (Farnell, Germany) and assembled into a female connector (Farnell, Germany). The lower half of the pins and all connecting wires were embedded in dental acrylate. To habituate the animal to the experimenter, regular handling was applied before the experiments.

Table 1

Stereotaxic coordinates of recording electrodes. Numbers in mm from bregma, according to Paxinos and Franklin (2004).

Regions	Coordinates (A/L/V)
OB ipsi	4.5/-0.8/-1.4
ACC contra	2/0.35/-1.1
ACC ipsi	2/-0.35/-1.1
S1 contra	−0.5/1.4/-0.7
S1 ipsi	−0.5/-1.4/-0.7
Ins contra	−1.1/3.7/-3.7
Ins ipsi	−1.1/-3.7/-3.7
VPL contra	−1.7/1.7/-3.35
VPL ipsi	−1.7/-1.7/-3.35
PAC ipsi	−2.7/-2/0
AMYG ipsi	−1.1/-2.3/4.0

Fig. 1

Schematic diagram of the experimental procedures.

A, Design of the experiment. Electrode implantation: Scheme of brain regions with electrodes chronically implanted ipsi- and/or contralateral to injection side: ACC: anterior cingulate cortex; S1: primary somatosensory cortex; Ins: insular cortex; VPL: ventral posterior-lateral thalamic nucleus; AMYG: central nucleus of amygdala; PAC: parietal cortex; OB: olfactory bulb; Two screws were implanted above the cerebellum as ground (GND) and reference (REF) electrodes. Habituation recordings: first local field potential (LFP) recordings in animal’s home cage for one hour (1 h) and four hours (4 h) the following day. Saline injection: One week later LFP recordings in home cage with saline (Sal) injections (10 μl) into the left hind paw. Two minutes (2 min) from the beginning of the recording (baseline: BL grey) and two minutes after Sal (blue) were used for data analysis. Capsaicin injection: one day later Capsaicin (Caps) (10 μl, 5 mg/ml) was administrated into the same paw the next day during the recording. Two minutes immediately after Caps (red) as well as the first two minutes one hour after Caps (green) injection were used for data analysis. B, C, D, Histology. Examples of Nissl stained slices showing the location of electrodes in ACC (B), S1 (C) and VPL (D). Red arrow indicates the position of the electrode tip. Scale bar: 1 mm; Right part of each panel indicates electrode positions indicated by red lines on schematic illustrations (Paxinos and Franklin, 2004).

Evaluating pain sensation

To evaluate pain sensation an electronic von Frey device (Bioseb, France) was used. For testing, each mouse was transferred into a vertical acrylic glass cylinder with a metal mesh on the base, so that the hind paw of the animal could be reached from below. The von Frey device was carefully placed on the plantar surface of the paw gradually increasing force until the animal retracted the paw. The test was performed three times for each hind paw alternating between the two sides. For analysis, the mean of these three tests was calculated.

Electrophysiological recordings

One week after surgery, first LFP recordings were performed in the animal’s home cage (Fig. 1). Electrophysiological signals were amplified (INTAN RHA2000 16 channels with a custom-added 3-d-accelerometer connected to three of the 16 input channels), band pass filtered (1–500 Hz, including the accelerometer signals), digitized (2.5 kHz) and stored on a hard disk for offline analysis. One week later, home cage recordings were repeated prior and immediately after saline (Sal) injections into the left hind paw. One day later 10 μl of 5 % capsaicin (Caps) solution was injected into the left hind paw (Fig. 1). The capsaicin solution (5 %) was prepared from 200 mg capsaicin powder solved in a mixture of 2 ml ethanol and 38 ml iso-tonic NaCl solution.

Data analysis

Electrophysiological Data was analyzed with built-in and custom-written MATLAB (Mathworks, USA) routines. Raw data was inspected visually and periods with artifacts were excluded from further analysis. For assessing pain states, recordings from the first two minutes directly after the respective injections (Sal, Caps) were compared to two minutes recorded under baseline (BL) conditions one hour before saline injection. Additionally, two minutes of recordings one hour after capsaicin injection (Caps (1 h)) were analyzed to test whether effects were reversible within the expected time. Periods with sleep were excluded from analysis.

Movement was assessed by three orthogonally mounted accelerometers on the head stage and was compared between the various control and pain conditions. Gravity-induced signals were removed by the band pass filter (see above, 1–500 Hz). Data was analyzed by dividing the integral of each of the three accelerometer channel signals by the epoch length of two seconds and then computing the mean of all three accelerometer channels over all epochs.

Respiration rate was calculated from the number of inspiratory and expiratory cycles per second recorded from nasal thermal sensors (see above) or from respiration-related oscillations in the olfactory bulb (Jessberger et al., 2016).

Power spectral density (PSD) was computed using the pwelch function (MATLAB, Signal Processing Toolbox) with a Hamming window of 4 s and 50 % overlap. Six frequency bands were defined as following: 1−4 Hz, 4−12 Hz, 12−30 Hz, 30−80 Hz, 80−120 Hz, 120−160 Hz. Band power for each frequency range was calculated by summing up all PSD values of the respective band. Power was normalized according to the following formula:

Coherence between brain regions was computed using the mscohere.m function (MATLAB Signal Processing Toolbox) for windows of two seconds with 50 % overlap for each of the six frequency bands. The group averages of interregional coherences for each frequency band were calculated and normalized to baseline with the same formula as the PSD values above

Cross-frequency coupling (CFC) was calculated using the method described by Tort et al. (Tort et al., 2010). In short, the raw signal was bandpass filtered for the phase (slow) and amplitude (fast) frequency band of interest. For each phase frequency band a phase time series was computed and for each amplitude frequency band an amplitude envelope time series was computed. The phase time series was then divided into phase bins of which each bin was related to the respective amplitude envelope value. Thus, a mean amplitude distribution for each phase bin was obtained. The divergence of this distribution from the uniform distribution was computed as a measurement of how strong the amplitude of the high frequency was modulated by the phase of the slow frequency (Modulation Index, ranging from 0 to 1). The modulation index was normalized to baseline with the same formula as the PSD values above.

Statistical analysis

Univariate tests were performed using a paired t-test, if not stated differently. In the case of non-Gaussian distributions, a Wilcoxon signed-rank test was applied. Due to the multiple testing problem, p-values were corrected by the procedure of Benjamini and Hochberg (Benjamini and Hochberg, 1995) to protect the false discovery rate. Respiration rate and movement were compared between groups using one-way ANOVA test for repeated measures corrected by the Tukey test of post-comparison. The von Frey hair test was compared between groups by two-way repeated-measures ANOVA with Bonferroni’s multiple comparison test. In the case of non-Gaussian distributions, we used the non-parametric Friedman test and Dunn's multiple comparisons test.

Due to the large number of hypotheses, finding a significant result by univariate testing was very unlikely. Therefore, a multivariate regression analysis was performed to detect potential differences when combining several parameters. To identify a set of parameters which allows differentiating between saline and capsaicin injection, multivariate logistic regression models were fitted. There were 32 observations (16 mice, each with saline and capsaicin injection) and 66 parameters for PSD (6 frequency bands, 11 regions), 11 parameters for CFC, and 330 parameters for coherence (6 frequency bands, 55 combinations of interregional coherence). Missing values of these parameters due to electrode malfunction were imputed by multiple imputation (van Buuren, 2018) based on a predictive mean matching approach (Rubin, 1986). Over all animals 14 recording channels had to be excluded resulting in 7.95 % of missing values for PSD and CFC analysis and 14.20 % of missing values for Coherence analysis. The region missing most often was S1i in three animals.

Since there are many more parameters than observations, a variable selection was appropriate for two reasons: Firstly, reducing the space of variables reduces noise in the data and thereby reduces the probability that conclusions are driven by chance. Secondly, selecting a suitable set of variables allows identifying those parameters that fit best to differentiate between saline and capsaicin injection. The elastic net (Zou and Hastie, 2005; Friedman et al., 2010) provides a well-suited method for variable selection in p > n situations, i.e., if there are more variables than data points. While in logistic regression the expectation of the outcome variable Y given the covariates X is modelled as , in the case of an elastic net this formula is restricted by the constraint . This implies that a sparser vector β will be chosen and therefore, fewer variables will be included in the model.

The parameters t and α are chosen by grid search and 8-fold cross-validation. Here, the dataset is split randomly into 8 equally sized folds. Then, for each combination of t and α out of the predefined grid, each fold is used exactly once as testing set and an elastic net is trained on the remaining 7 folds in order to predict the missing one. The mean accuracy over the resulting 8 predictions is applied as metric to evaluate which combination of t and α allows the best prediction. Finally, the parameters with a non-zero value of β are included in the model.

For each analysis, the selected variables (PSD, CFC, and coherence) were fitted to a logistic regression model. The outcome variable was the type of injection (saline or capsaicin, i.e., control or pain) and the selected variables were included as covariates. We adjusted for the fact that each mouse was measured twice by incorporating the random effect “mouse” in the model. Receiver operating characteristics (ROC) is a commonly used tool for evaluating the performance of classification models. For each possible threshold to classify between two outcomes the respective sensitivities and specificities are calculated and plotted with the sensitivity as ordinate and the specificity as abscissae. The area under the curve then represents the performance of the classification with 1 being a perfect classification and 0.5 corresponding to chance level. ROC-curves were plotted and the corresponding AUC value, as well as the prediction accuracy were computed with 95 % confidence intervals. All confidence intervals (CIs) are computed by means of the formula by Wilson (Wilson, 1927). To test whether the AUC differed significantly from chance (AUC = 0.5) a p-value was derived by a Wilcoxon-Mann-Whitney U test (DeLong et al., 1988). To avoid overfitting and to show a realistic estimate of these curves, they were also created by 8-fold cross-validation and therefore beta values were computed as a mean over all cross-validations. Univariate tests for significance were performed in MATLAB. All further calculations were done with the statistical software R (R Core Team, 2018).

Histology

After completion of the experiments the animals were anesthetized with isoflurane (see above) and direct anodal current (20 μA) was passed for 10−15 s to one wire of each of the depth electrodes for marking the localization of the electrode positions. Three days later, the animals were deeply anesthetized with phentobarbital, perfused with polychlorinated biphenyl (PCB) followed by paraformaldehyde (PFA). The brain was carefully removed and stored in PFA solution at 4 degrees Celsius. The fixated brain was cut in slices of 50 μm using a Vibratom (Leica VT 1200S). Slices were dried and stained with Cresyl Violet (Nissl staining) before confirming the correct electrode position by microscopy (see Fig. 1).

Results

Data are based on analysis of two-minute periods of local field potentials (LFP) recorded from selected brain regions of freely moving mice during four conditions: first, prior to any manipulation (baseline, BL); second, immediately after injection of saline into the left hind paw (Sal); third, immediately after injection of capsaicin (Caps); and fourth, one hour after injection of capsaicin (Caps (1 h)). See Fig. 1, Table 1 and Methods for further details.

Fig. 2 shows examples of raw signals in olfactory bulb (OB), anterior cingulate cortex (ACC), insular cortex (Ins), primary somatosensory cortex (S1), ventral posterolateral thalamic nucleus (VPL), central nucleus of the amygdala (AMYG) and the parietal cortex (PAC) during the four conditions defined above. Note that recordings were suited to include the respiration-entrained rhythm (RR) which overlaps with theta oscillations (Nguyen Chi et al., 2016) and covers frequencies up to ∼8 Hz in the present series of experiments. RR is most pronounced in the olfactory bulb but also clearly visible in ACC, Ins and VPL.

Fig. 2

Raw traces of local field potentials for selected brain areas and experimental conditions.

Brain areas: OB: olfactory bulb; ACC: anterior cingulate cortex; Ins: Insular cortex; S1: primary somatosensory cortex; VPL: ventral posterior-lateral thalamic nucleus; AMYG: central nucleus of amygdala; PAC: parietal cortex. ACC, Ins, S1 are VPL are contralateral to the injection side, OB, AMYG and PAC are from the ipsilateral hemisphere. First column: baseline (BL) recording during first two minutes one to two hours prior to saline injection (Sal); second column: first two minutes after Sal injection; third column: first two minutes after capsaicin (Caps) injection; fourth column: first two minutes one hour after Caps injection (Caps (1 h)). Note the presence of low frequency oscillations in all brain regions which follow the OB. Recordings were taken from the hemisphere contralateral to the injection side with exception of OB, AMYG and PAC which were recorded ipsilaterally.

To exclude respiration rate and movement activity as confounding factors we compared mean respiration rates (Fig. 3A) and movement activity, based on accelerometer integrals (Fig. 3B) between the four analyzed conditions. Respiration rates differed between the four conditions (one-way ANOVA test; F (3,45) = 4.065; p < 0.05; n = 16). Multiple comparisons revealed that there was a significant decrease of respiration rate between baseline (BL) and Caps (from mean values of 8.5 Hz to 7.5 Hz, p < 0.05; n = 16; Tukey's multiple comparisons test, F (3,45) = 4.065, p < 0.05, n = 16)) but there was no difference in any other pairwise comparison, most importantly not between control (Sal) and Caps. No differences in integrated motor activity were found between the four conditions (one-way ANOVA test; F (3,45) = 0.2427; p > 0.5; n = 16; Fig. 3B).

Fig. 3

Enhanced hyperalgesia after capsaicin compared to saline but no difference in respiration rate or movement.

A: Mean Respiration rate with SEM during baseline (BL), the first two minutes after saline (Sal) or capsaicin (Caps) injections as well as one hour after capsaicin injection (Caps (1 h)). There was a significant decrease of respiration rate between baseline (BL) and Caps (from mean values of 8.5 Hz to 7.5 Hz). Note, there is no difference between Caps and the control Sal. One-way ANOVA with Tukey's multiple comparisons test, F (3,45) = 4.065; *p < 0.05; n = 16. B: Mean movement based on accelerometer integrals. There is no difference between the conditions. One-way ANOVA; F (3,45) = 0.2427; p > 0.5; n = 16. Abbreviations as in A. C: Enhanced pain sensation two minutes (2 min; *** p < 0.001) and one hour (1 h; ** p < 0.01) after Caps injection based on the withdrawal threshold during the von Frey hair test (Two-way repeated-measures ANOVA with Bonferroni’s multiple comparisons test, n = 6.).

Visual inspection revealed that mice frequently licked the left hind paw during the first two minutes after capsaicin injection, indicating ongoing pain. There was no licking of the paw after saline injection. In six animals pain sensation was separately verified using the electronic von Frey test which measures the threshold of paw withdrawal upon mechanical stimulation (for further details see methods). Fig. 3C shows the significant difference between the control (saline) group and the capsaicin group (two-way ANOVA; F(1,15) = 40,71; p < 0.0001; n = 6). Mean withdrawal threshold significantly decreased at two minutes (2 min) and one hour (1 h) after Capsaicin injection (compared to Sal injection; p < 0.001 at 2 min and p < 0.01 at 1 h; n = 6; Bonferroni’s multiple comparison test). Thus, there were clear indications for tonic pain following capsaicin injection. Hyperalgesic responses lasted for at least one hour after the nociceptive cue.

Next, we analyzed the power spectral density in all frequency bands and all measured brain regions during two-minute periods in the four states chosen for analysis. Power was normalized according to the formula shown in Methods. No power changes were found in acute pain (Caps: the first two minutes after capsaicin injection) compared to saline (Sal) injection when correcting for multiple comparisons (without correction, the ipsilateral ACC showed an apparent power decrease in the 80−120 Hz band; p < 0.05)

In a next step we analyzed cross-frequency coupling (CFC) between the respiration-related low frequency rhythm (RR) and sub-bands of gamma. We measured the modulation between low-frequency phase (1–15 Hz, abscissae) and high-frequency amplitude (20–200 Hz, ordinate). Without correction for multiple comparisons, cross-frequency coupling between 1−12 Hz and 80−120 Hz was enhanced with p < 0.05 in four regions (ACCc, ACCi, S1i and PACi). Fig. 4A exemplary shows CFC in the ipsilateral anterior cingulate cortex (ACC ipsi). Coupling between slow and fast oscillations increased markedly after capsaicin injection and declined within 1 h (Fig. 4B). At the same time, absolute values of power for gamma rather decreased (80−120 Hz, Fig. 4C) and slow-frequency (1−12 Hz) oscillations (Fig. 4D) did not change between Sal and Caps. This suggests that the observed increase in phase-amplitude coupling was not secondary to a mere increase in the respective power bands. After correcting for multiple testing, however, none of these parameters remained significantly different from saline injection.

Fig. 4

Cross-frequency coupling in ipsilateral anterior cingulate cortex after capsaicin-induced pain.

A: Cross-frequency coupling (CFC) in ipsilateral anterior cingulate cortex (ACC ipsi) two minutes after injections of saline (Sal) or capsaicin (Caps) and one hour after Capsaicin (Caps 1 h). Averaged heat maps (n = 14) are shown. Ordinates: Amplitudes for frequencies from 20 to 200 Hz; Abscissae: Phases for frequencies between 1 and 15 Hz. The color depicts the modulation index (MI) between oscillations of both frequencies (see Methods) (Tort et al., 2010). Superimposed white solid lines indicate the power spectral density (PSD) of local field potentials of the respective brain regions and, yellow dashed lines correspond to the LFP in OB corresponding to respiration rate (see text) (Yanovsky et al., 2014). Note increased CFC in Caps B: Averaged modulation index (MI) between phases of low frequencies (1–12 Hz) and amplitudes of gamma (80–120) as shown in A differ between Sal (blue) and Caps (red) but recover one hour later (green). MI values were normalized to baseline prior to averaging. Means with SEM. Multiple t-test, *p < 0.05 (without multiple corrections); n = 14. C: Mean gamma (80−120 Hz) band power decreased between Sal and Caps. D: Mean power of low frequencies (1−12 Hz) did not change between Sal and Caps (For C and D, Means with SEM. Multiple t-test, *p < 0.05, **p < 0.01 (before multiple corrections); n = 14).

Pain processing occurs in a widely distributed network. We therefore looked at coherence of rhythmic activity, making use of our parallel recordings from multiple brain areas and frequency bands. Coherence was calculated between all eleven recorded brain regions for all six frequency ranges after Sal and Caps injection, respectively. Several differences between Sal and Caps conditions were found, mostly at frequencies up to 30 Hz. Fig. 5A shows color-coded maps of coherence coefficients for Sal (upper row) and Caps (second row) in three frequency ranges: 1−4 Hz (left column), 4−12 Hz (middle) and 12−30 Hz (right). Apparent differences between Caps and Sal with p < 0.05 are shown in Fig. 5B in a visualization of their topographic relation. Most interregional coherence differences were found in the frequency range between 4−12 Hz. Coherence was increased under Caps conditions between OBi and VPLc / VPLi; between ACCi and ACCc / S1i / VPLc / VPLi; between ACCc and Insi / AMYGi / VPLc / VPLi; between S1c and Insi / AMYGi / VPLc / VPLi; between S1i and Insi / AMYGi; between OBi and VPLc / VPLi. In the frequency range 12−30 Hz coherence increased between PACi and OBi / VPLi; and between S1c and Insc. The only decrease in coherence after capsaicin injection (compared to Sal) was found between S1c and Insc. In all three higher frequency ranges there were only few changes in interregional coherence: Following capsaicin injection coherence increased in the 30−80 Hz range between OBi and AMYGi / VPLi; in the 120−160 Hz range coherence increased between OBi and Insi. None of the apparent pairwise changes, however, remained significant after correcting for multiple comparisons.

Fig. 5

Coherence between certain brain regions during acute pain mainly in frequencies up to 30 Hz.

A: Color coded maps of coherences between eleven brain regions: ipsilateral olfactory bulb (OBi: number code 1), ipsi- and contralateral anterior cingulate cortex (ACCi: 3 and ACCc: 2), ipsi- and contralateral primary somatosensory cortex (S1i: 5 and S1c: 4), ipsi- and contralateral insular cortex (Insi: 7 and Insc: 6), ipsilateral ventral posterolateral thalamic nucleus (VPLi: 10 and VPLc: 9) and ipsilateral parietal cortex (PACi: 11). The left column are coherences for frequencies from 1 to 4 Hz, the middle column for frequencies from 4 to 12 Hz and the right column for frequencies from 12 to 30 Hz. The first row shows the changes of coherences after saline injection compared to baseline (Sal), the second row illustrates the changes after capsaicin injection compared to baseline (Caps) and the third row displays the difference between Sal and 1 h after Capsaicin injection(Caps 1 h) injection. B: Schematic summary of differences with p < 0.05 in interregional coherences for the three frequency bands for Caps compared to Sal. C: Schematic summary of differences with p < 0.05 in interregional coherences for the three frequency bands for Caps 1 h compared to Sal.

Fig. 5C shows apparent pairwise changes of coherence one hour after capsaicin injection compared to saline injection with p < 0.05. In the 1−4 Hz frequency band coherence decreased between AMYGi and ACCc / PACi; between OBi and PACi; between S1i and ACCc. An increase in coherence was found between S1i and Insi. In the 4−12 Hz band the only remaining change was an increase in coherence between ACCc and VPLi. In the 12−30 Hz frequency band coherence increased between OBi and Insi / VPLi / AMYGi / PACi; between Insi and VPLi / VPLc / AMYGi / PACi; between AMYGi and VPLi / PACi; between S1c and VPLc. The 30−80 Hz band showed a similar pattern with increased coherence between OBi and Insi / VPLi / PACi; between Insi and AMYGi / S1i / PACi; between VPLi and AMYGi / S1i; between S1c and ACCc / AMYGi. In the 80−120 Hz band coherence increased between ACCc and AMYGi. In the 120–160 HZ band coherence increased between OBi and Insi; between ACCc and PACi / VPLc / AMYGi. Again, none of these apparent changes remained significant after correcting for multiple testing.

In summary, oscillation coherence seems to increase between multiple brain regions under tonic pain conditions, most markedly in the frequency range of 4–12 Hz. Fig. 6B shows, as an example, raw values of coherence between OBi and VPLi for the frequency range of 4–12 Hz. Most of these apparent changes seem to recover one hour after capsaicin injection, but at this time there were remarkable increases of coherence in the 12−30 Hz and 30−80 Hz frequency range, mostly involving OBi and other ipsilateral regions. The high number of pairwise combinations, however, limits the power of simple pairwise comparisons which did not remain statistically significant following correction for multiple comparisons.

Fig. 6

Absolute coherence values.

A: Histogram showing the frequency of occurrence of group mean coherence values for all frequency bands, combinations of regions, and both saline and capsaicin injection. The minimum value of coherence from all pairwise tests was 0.031, the maximum value was 0.836. B: Example of absolute group mean coherence values between OBi and VPLi for saline (black) and capsaicin injection (red) in the theta frequency band (4−12 Hz). Shaded regions represent standard error of the mean.

Changes of distributed brain activity may be best detected as complex patterns, rather than as multiple differences in individual parameters. In order to uncover such characteristic patterns under pain we next employed elastic net analysis as a variable reduction method (Zou and Hastie, 2005). With this method the most influential parameters for classifying pain versus non-pain stats can be defined without referring to the significance of individual pairwise comparisons (see Methods). As a result, out of sixty-six potential variables from power analysis, fourteen values were selected. For these parameters, the area under the curve (AUC) in the receiver operating characteristic (ROC) analysis was 0.650 (95 % CI [0.477, 0.791]) (Fig. 7A), not significantly different from 0.5 (chance level). Likewise, the use of seven (out of eleven) selected CFC variables (Fig. 7B) did not permit a classification between saline and capsaicin injections, despite an AUC of 0.711 (95 % CI [0.538, 0.838]). The nine (out of 330) coherence values alone (Fig. 7C), however, did result in a classification that was better than chance (AUC 0.785 (95 % CI [0.617, 0.893])), selected variables were S1i-VPLi (4−12 Hz, β: 80.1), OBi-ACCc (4−12 Hz, β: 114.3), OBi-VPLc (4−12 Hz, β: 124.4), OBi-VPLi (4−12 Hz, β: 191.5), OBi-VPLi (12−30 Hz, β: 46.9), OBi-PACi (12−30 Hz, β: 62.7), S1c-Insc (12−30 Hz, β: 195.0), OBi-Insi (30−80 Hz, β: 15.4), Insi-PACi (30−80 Hz, β: 138.8). Next, ten variables were selected from a combination of all three types of data, i.e. power, CFC and coherence (Fig. 7D). This selection resulted in a classification with an AUC of 0.926 (95 % CI [0.783, 0.977]), significantly different from chance. Selected variables were Coherence of S1i-VPLi (4−12 Hz, β: 66.4), OBi-ACCc (4−12 Hz, β: 121.1), OBi-VPLc (4−12 Hz, β: 144.5), OBi-VPLi (4−12 Hz, β: 184.4), OBi-VPLi (12−30 Hz, β: 0.4), OBi-PACi (12−30 Hz, β: 51.1), S1c-Insc (12−30 Hz, β: 160.3), Insi-PACi (30−80 Hz, β: 99.1) and PSD of AMYGi (120−160 Hz, β: 35.1). Note that no CFC variable contributed to the final selection. A decreased, but still better than chance, performance was achieved when only using parameters with a frequency below 30 Hz. In this case, five (out of 198) parameters were selected with the algorithm (Fig. 7E, AUC = 0.738 (95 % CI [0.567, 0.859]), significantly different from chance). Selected parameters were coherence between S1i-VPLi (4−12 Hz, β: 9), OBi-VPLc (4−12 Hz, β: 2.0), OBi-VPLi (4−12 Hz, β: 9.6), OBi-VPLi (12−30 Hz, β: 12.0) and PSD of S1c (12−30 Hz, β: 4.9). To test for recovery, we also tried to classify between the states at saline and one hour after capsaicin injection, respectively. Using the same methods as above we could not find any classification with results better than chance, except for CFC (AUC = 0.719 (95 % CI [0.546, 0.844], selected regions: S1c, ACCi, VPLc, VPLi, Insc, Insi, PACi).

Fig. 7

Classification analysis after variable selection with elastic net model predicts pain with parameters up to 30 Hz.

The Elastic net model was used to select variables for each analysis respectively and then fit a logistic regression with these variables to classify between saline and capsaicin injection. The classification performance was evaluated by computing the area under the curve (AUC, between the black and gray lines) of a receiver operating characteristic (ROC). A: For power spectral density (PSD) 14 out of 66 variables were selected resulting in an AUC of 0.650. B: For cross-frequency coupling (CFC) 7 out of 11 variables were selected resulting in an AUC of 0.711. C: For interregional coherences 9 out of 360 variables were selected resulting in an AUC of 0.785. D: For all three analyses (PSD, CFC and Coherence) together 9 out of 407 variables were selected resulting in an AUC of 0.925. E: Using only frequencies below 30 Hz for PSD and Coherences, 3 out of 198 variables selected resulting in an AUC of 0.738. F: Comparison between each ROC curve for measurements taken immediately after injections and at 1 h after capsaicin injection, respectively, using the DeLong test. All p-values were corrected for multiple testing with the procedure by Benjamini and Hochberg.

We also tested for unspecific effects of sex, respiration rate or movement (‘clinical’). These parameters were included in all classification analyses but were not selected as variables. Using sex, respiration rate and movement alone, no classification was possible (AUC = 0.566 (95 % CI [0.397, 0.722]). Moreover, when they were included in the best model the AUC decreased from 0.926 to 0.801. A Chi-square test using sex and injection type was not significant. To statistically compare the ROC curves, we used the DeLong test which showed a p-value below 0.05 for clinical vs. all. This indicates that using all electrophysiological parameters showed a better classification than using the clinical parameters (sex, respiration rate or movement, Fig. 7F). For the ROC curves one hour after injection there was no p-value below 0.05. Again, all p-values were adjusted by the procedure by Benjamini and Hochberg to correct for multiple testing.

Discussion

We report capsaicin-induced changes of electrophysiological signals at the network level in a wide range of brain regions in the mouse. We found several changes of electrical activity following injection of capsaicin as compared to saline: Cross-frequency coupling between low (1−12 Hz) and high (80−120 Hz) frequency oscillations was increased in several brain regions, as was the power within a sub-band of gamma (80–120) in ipsilateral ACC. Similarly, interregional coherence between different brain areas was increased, mainly in the theta frequency range (4−12 Hz). However, due to the large number of parameters being compared, none of these changes was significant after controlling for multiple testing. To overcome this problem we performed a variable selection with an elastic net model (Zou and Hastie 2005). With as few as ten selected variables a classification of pain versus non-pain with an AUC of 0.926 could be achieved. When analysis was confined to frequencies below 30 Hz, the elastic net approach still found a significant set of five parameters, resulting in an AUC of 0.738. While the present study is restrained to data from mice within a specific experimental paradigm, similar approaches may be helpful in the search for useful classifiers within human routine EEG frequencies.

In contrast to previous reports from rats and humans (Leblanc et al., 2014, 2015; LeBlanc et al., 2016; Nickel et al., 2017) we did not find an increase of power in theta frequency or other frequency ranges during acute pain compared to saline injection. This might reflect differences in species (rats versus mice) and methods of pain induction (laser-evoked pain versus capsaicin). In any case, the lack of clear power changes within standard frequency bands led us to probe higher-order phenomena such as cross-frequency coupling and inter-regional coherence. Within these parameters, we found a capsaicin-specific modulation of gamma oscillations in the 80–120 Hz band by slow (1–12 Hz) oscillations in the ipsilateral and contralateral anterior cingulate cortex (ACC) as well as ipsilateral somatosensory cortex (S1i) and parietal cortex (PACi). Within this frequency range, several types of network activity overlap including sleep-related slow oscillations, theta oscillations and respiration-related network oscillations (Yanovsky et al., 2014; Tort et al., 2018). Our data were recorded from awake mice exhibiting theta oscillations, such that sleep-related slow oscillations are unlikely to occur. Theta and respiration-related oscillations, however, occur simultaneously and have overlapping frequencies. Interestingly, the gamma frequency band (80−120 Hz) with altered coupling to slow frequencies is selectively modulated by RR, in contrast to theta oscillations (Zhong et al., 2017). It is thus likely that the pain-related modulation of CFC is linked to respiration.

An ongoing discussion regarding slow oscillations in the rodent brain is whether these are really generated by neuronal mechanisms at the location of measurement or passively transferred over considerable distances by mere biophysical processes (Kajikawa and Schroeder, 2011). The monopolar electrodes used in our present recordings do not allow distinction between both mechanisms. While our classifier would not lose validity even if oscillations were volume-conducted to the respective sites, we postulate that slow oscillations are indeed generated or accompanied by clear neurophysiological mechanisms at the sites of measurements. We have previously shown that respiration-related oscillations i) have defined current sinks and sources in the mouse hippocampus (Yanovsky et al., 2014), ii) selectively engage a specific gamma subband in CFC (Zhong et al., 2017), iii) rhythmically entrain unit discharges (i.e. neuronal spikes) in hippocampal, parietal and prefrontal cortex (Yanovsky et al., 2014; Nguyen et al., 2016; Zhong et al., 2017; Tort et al., 2018). Similar findings did also apply to theta oscillations, irrespective of the fact that these might indeed be volume-conducted from the hippocampus to adjacent cortical areas (Sirota et al., 2003). From these data we can state that slow rhythmic network activity as measured in the present study is at least not exclusively mediated by volume conduction, and that wide-ranging slow rhythms have direct effects on neuronal discharges in at least some of the regions under study.

At the same time, our present findings underline the importance of the anterior cingulate cortex (ACC) in pain processing (Hutchison et al., 1999; Lieberman and Eisenberger, 2015; Bliss et al., 2016; Meda et al., 2019). Interestingly we found increased CFC in both, the contralateral and ipsilateral ACC. As most pain pathways are crossing sides in the spinal cord this finding points towards a higher-order process in the anterior cingulum, in line with its known role in attention and emotional processing of pain (Rainville, 1997; Shackman et al., 2011). This is supported by observations showing that some neurons in the cingulum respond to ipsilateral and contralateral painful stimuli (Hutchison et al., 1999) and that neuropathic pain causes metabolic changes primarily in the ipsilateral cingulate cortex, amongst other structures (Shen et al., 2019).

Regarding interregional coherence, changes were even more pronounced between ipsilateral regions. Coherence is thought to reflect communication between distant brain regions (Fries, 2015) and thus may play an important role in forming an integrated perception of acute pain. Interestingly, increases of coherence were mostly restricted to frequencies below 30 Hz. This underlines the important role of slow oscillations for long-range coordination of network activity (Reinhart and Nguyen, 2019), and – if applicable to humans – it facilitates the search for clinically usable diagnostic parameters. Indeed, previous studies on acute pain in humans revealed a suppression of frequencies within the alpha range (Nir et al., 2012; Shao et al., 2012; Li et al., 2016b). Yet, such alterations are a common effect associated with many mental processes, while pain-related changes in the gamma frequency range (Peng et al., 2014; Schulz et al., 2015; Li et al., 2016a) might be more specific. Likewise, higher-order parameters of network activity (CFC, coherence) reflect complex dynamic interactions within the pain matrix and may, hence, provide more specific signatures of pain (Kucyi and Davis, 2015).

In any case, it remains challenging to discriminate physiological markers of genuine pain from saliency (Legrain et al., 2011). This does also apply to the saline controls in our present experiments. We note, however, that capsaicin induced a clear hyperalgesic reaction (lower threshold for paw withdrawal upon mechanical stimulation) together with typical behavioral signs (paw licking). Such behaviors were absent after saline injection and are, thus, indicative of capsaicin-induced pain. In line with these behavioral data, our electrophysiological classifiers were able to distinguish between capsaicin and saline at early times after the injection, but not after one hour (see Fig. 7). The underlying electrophysiological data (power, CFC, coherence) were normalized to values measured in the same environment before any injection, and hence indicate specific changes in reaction to one substance (capsaicin) but not the other (saline). Together, it is thus very likely that the observed changes in the early phase following capsaicin-injection were induced by pain, rather than by unspecific changes following the injection.

In order to develop a practicable method for distinguishing pain from control states we used an elastic net model for variable reduction and selection of relevant parameters. Such approaches are necessary to avoid overfitting of a limited data set and to define small sets of parameters which are easily applicable in clinical settings. Our analysis did indeed reveal informative sets of five to nine parameters allowing a distinction with AUC values up to > 0.9. Interestingly, this analysis used almost exclusively coherence parameters, except one PSD value from AMYG ipsi (120−160 Hz). Thus, inter-regional coherence of network oscillations is a key feature of network activity following acute pain induction. Finally, when using only oscillations at frequencies below 30 Hz our algorithm could still identify a reasonable classification (AUC 0.738) which might be useful for similar approaches in EEG recordings of limited bandwidth, e.g. in humans.

In conclusion, it seems promising to further investigate pain-specific signatures of acute pain using electrophysiological measures of brain activity in freely moving mice. In our paradigm, changes were mostly present in interregional coherence underlining that pain is processed in a distributed, yet co-active brain-wide matrix. The underlying mechanisms and the significance of each parameter for pain processing remain to be elucidated (Mouraux and Iannetti, 2018). It is also not clear whether the observed alterations do reflect the genuine experience of pain or result from secondary changes in brain activity during the course of the experiment. Nevertheless, our data show that a combined set of (mostly higher-order) parameters of network activity is predictive of the state of capsaicin injection in mice. This finding may be useful in the ongoing search for a neurophysiological signature of pain at the network level.

Ethical standards

We comply with all ethical standards regarding animal experiments as indicated in the Methods section.

Conflicts of interest

The authors have no conflict of interest to declare.

CRediT authorship contribution statement

Simon Ponsel: Investigation, Writing - original draft, Writing - review & editing, Visualization. Jiaojiao Zhang: Investigation, Writing - original draft, Writing - review & editing, Visualization. Maximilian Pilz: Formal analysis, Validation, Visualization. Yevgenij Yanovsky: Conceptualization, Methodology, Supervision. Jurij Brankačk: Methodology, Supervision, Writing - original draft, Writing - review & editing. Andreas Draguhn: Conceptualization, Supervision, Project administration, Resources, Funding acquisition.