Characteristics of Evoked Potential Multiple EEG Recordings in Patients with Chronic Pain by Means of Parallel Factor Analysis

Abstract

This paper presents an alternative method, called as parallel factor analysis (PARAFAC) with a continuous wavelet transform, to analyze of brain activity in patients with chronic pain in the time-frequency-channel domain and quantifies differences between chronic pain patients and controls in these domains. The event related multiple EEG recordings of the chronic pain patients and non-pain controls with somatosensory stimuli (pain, random pain, touch, random touch) are analyzed. Multiple linear regression (MLR) is applied to describe the effects of aging on the frequency response differences between patients and controls. The results show that the somatosensory cortical responses occurred around 250 ms in both groups. In the frequency domain, the neural response frequency in the pain group (around 4 Hz) was less than that in the control group (around 5.5 Hz) under the somatosensory stimuli. In the channel domain, cortical activation was predominant in the frontal region for the chronic pain group and in the central region for controls. The indices of active ratios were statistical significant between the two groups in the frontal and central regions. These findings demonstrate that the PARAFAC is an interesting method to understanding the pathophysiological characteristics of chronic pain.

1. Introduction

Chronic pain is a complex disease characterized by pain persisting after damage or pathology has healed. Effective treatment of chronic pain is hampered by an incomplete understanding of the pathophysiological changes that occur in the nervous system of chronic pain sufferers. The electroencephalogram (EEG) records the electrical activity from the scalp produced by the firing of neurons within the cerebral cortex [1] and has been widely used to analyze neural activity in chronic pain subjects [2, 3].

To investigate the physiological basis of chronic pain, event-related potentials (ERPs) have been used to explore pain-related modulation of the latency, location, amplitude, and frequency of evoked EEG responses to sensory stimulation [4–6]. Previous studies have shown that frequency and time domain characteristics of EEG recordings from certain brain regions are altered with chronic pain [7–9]. Traditionally, the ERP components are analyzed at specific cortical locations, that is, the vertex and frontotemporal region (e.g., [10, 11]). With this approach, however, the information provided is limited to the particular region under investigation and neglects the importance of wider cortical regions in information processing [12].

Multiple EEG recordings can be used to characterise the electrical activity across the whole cortex. Traditional methods of PCA and ICA analysis process the multiple EEG signals in two-way domains such as time channel and frequency channel [13]. To effectively characterise multiple EEG signals, development of an analysis method that captures time-frequency-channel information is required.

In this study, we have investigated differences in multiple ERPs between chronic pain patients and pain-free individuals. To characterise the EEG signals in the time-frequency-channel domain, a parallel factor analysis (PARAFAC) method with wavelet transforms was developed to decompose the multiple EEG recordings. The PARAFAC method has been successfully employed to detect abnormal EEG activity in neurological diseases such as epilepsy and Alzheimer’s disease [14]. Herein, this novel method is employed to analyze the EEG in chronic pain subjects. The results show differences in cortical evoked activity between chronic pain and pain-free individuals and indicate that the PARAFAC is an effective method for extracting the characteristics of multiple EEG recordings in the time-frequency-channel domain.

2. Materials and Methods

2.1. Subjects

Subjects were 13 chronic pain patients recruited through the Waikato Hospital Pain Clinic and 13 pain-free volunteers. A range of conditions were represented in the patient group, including chronic lower back pain, neck pain, abdominal pain, and throat pain. None of the subjects in either group had a history of other neurological disease or head injury. Ethical approval was obtained from the Waikato Ethics Committee, and all subjects signed written informed consent.

2.2. Sensory Stimulations

The stimuli consisted of brief (10 ms), repetitive (at least 120) electrical shocks delivered to the dominant index finger. The electrodes were positioned on the dorsal aspect of the distal interphalangeal joint and a fine wire with a soldered tip, applied to the pulp of the fingertip. The stimulus intensity was recorded as the percent maximum voltage and rated by each subject on a 1–10 analogue scale. Two intensities of electric shock were tailored to each subject, one that was easily felt but not painful (from here on referred to as the “touch” stimulus) and one that was rated as “moderately painful” (the “pain” stimulus). The pain stimulus was felt as a sharp pricking sensation, predominantly under the wire electrode on the finger pulp. The shocks were given in three sequences at a constant frequency of 1 every 1.5 s as follows: (1) 120 sequential touch shocks (“touch”); (2) 120 sequential painful shocks (“pain”); (3) a random sequence of 300 touch and painful shocks at a 4 : 1 ratio. For the randomized protocol, the touch and pain stimuli were analyzed separately (“random touch” and “random pain”, resp.). Stimulus intensity and delivery were controlled by MatLab software (Matlab 6.0 Mathworks, Natick, MA, USA) running on a laptop computer that interfaced directly with the stimulus generator.

2.3. EEG Recording and Experimental Protocol

The subjects were comfortably seated, and the stimulating and EEG recording electrodes were attached. The latter consisted of a 28-channel bipolar montage configured in accordance with the international 10 : 20 system. The electrodes were Ag/AgCl sintered ring electrodes (Falk Minow, Herrsching, Germany) (1 cm outer diameter) that fastened securely to plastic loops imbedded in a prefabricated scalp cap (Easycap, Falk Minow, Herrsching, Germany). One of two cap sizes was chosen to give the correct positioning of the electrodes on the head relative to the nasion and inion, in accordance with the international 10 : 20 system. The centres of the electrodes were filled with an electrolyte gel, and attention was given to ensure the gel made contact with the scalp. Two reference electrodes were positioned behind each ear. The EEG electrodes were connected to two 16-channel biosignal amplifiers (Guger Technologies, Herbersteinstrasse, Austria) and digitised (Gdaqsys, Guger Technologies, Herbersteinstrasse, Austria) to computer at 100 Hz for continuous display and later offline analysis. The amplifiers were powered using mains-charged battery packs. One of the spare channels on the amplifier was used as an event marker from the electrical stimulus generator. Application of the electrodes took approximately 1 hour. The quality of the EEG was assessed by visual inspection and corrective measures taken to improve the quality of “noisy” channels. This usually involved checking the contact of the electrolyte gel between the electrode and the scalp. Time restraints, particularly with the requirement for patients to be seated for up to two hours to complete the study, meant it was not practicable to check and monitor individual channel impedances.

During delivery of the stimulation sequences, the subjects were instructed to keep their eyes closed, refrain from talking, and relax as much as possible. The subjects were not specifically instructed to either attend to or ignore the stimuli. The subjects could stop the stimulation at any time by pressing a button. The three sequences took approximately 30 minutes to complete.

2.4. Data Analysis

2.4.1. Preprocessing

The EEG was preprocessed by a band-pass filter and further analyzed using EEGLAB [15] software in MatLab. The raw EEG was 1 to 50 Hz band-pass filtered. Each trace was visually inspected, and data predominated by electrical noise was discarded.

2.4.2. Wavelet Transforms

Wavelet transform was used to transform a single-channel EEG signal into a time-frequency map. In this study, the continuous wavelet transform (CWT) was applied, and the Morlet wavelet was employed [16], and it is𝜓₀(𝑡)=𝜋^−1/4𝑒^𝑖𝑤𝑡𝑒^−1/2𝑡2, (1) where 𝑤 is the wavelet central angle frequency, often 𝑤≥6, which is an optimal value to adjust the time-frequency resolution [17]. In this study, 𝑤=6 was applied. Then, a family of wavelets can be generated: 𝜓_s(𝑡)=(1/√2)𝜓₀(𝑡/𝑠), s∈(0,+∞), and 𝑠 is called a scale. The CWT at scale 𝑠 and time 𝑡 of a signal 𝑥(𝑡) is defined as𝑊()=1𝑠,𝜏√𝑠𝑥(𝑡)𝜓_𝑠∗𝑡−𝜏𝑠𝑑𝑡, (2) where 𝜓_𝑠 is a parent wavelet function, and ∗ denotes complex conjugation. By adjusting the scale 𝑠 and the translation 𝜏, a series of different frequency resolutions in the signal can be projected on the two-dimension space (scale 𝑠 and translation 𝜏). The factor √𝑠 normalizes energy across the different scales.

2.4.3. PARAFAC

After wavelet transformation of all EEG channels, a three-way tensor 𝑋(𝑡, 𝑓, 𝑐) (time-frequency-channel), giving the energy at time 𝑡, frequency 𝑓, and channel 𝑐, was obtained. To decompose the three-way tensor into time, frequency, and channel modes, the PARAFAC method was applied, and a linear combination of the three-way tensor was obtained by means of the alternating least squares (ALS) algorithm [14, 18]. The PARAFAC model is defined as𝑋_{𝑇×𝐹×𝐶}=^N_𝑛=1𝑎_𝑛∘𝑏_𝑛∘𝑐_𝑛+𝐸, (3) where 𝑁 is the number of signal factors of 𝑋_{𝑇×𝐹×𝐶}, and 𝑎_𝑛, 𝑏_𝑛, and 𝑐_𝑛 indicate the 𝑛th column of the loading matrices 𝐴∈𝑅^𝑇×𝑁, 𝐵∈𝑅^𝐹×𝑁, and 𝐶∈𝑅^𝐶×𝑁, respectively. 𝐴, 𝐵, and 𝐶 represent the time, frequency, and channel modes and provide information on the interactions between modes. 𝐸∈𝑅^{𝑇×𝐹×𝐶} is the residual information in the decomposition. The operator ∘ represents the outer product of two vectors. Illustration of a 2-factor PARAFAC model on a three-way dataset is shown in Figure 1(A).

Figure 1 

The PARAFAC model and factors extracted by the PARAFAC model from one case. (A) PARAFAC modeling of a three-way tensor. Each component (𝑅=2) is the outer product of 𝑎, 𝑏, and 𝑐 of rank-1, and E is a residual tensor. (B) Flowchart describing the procedure of the epoch decomposition process by the PARAFAC method. (C) Whole brain cortex is clustered into 5 zones (1 = frontal, 2 = central, 3 = occipital, 4 = left temporal, and 5 = right temporal).

In the PARAFAC method, determination of the number of factors is a key issue. There are several methods to determine the number of factors, including the visual appearance of loadings, the residual analysis, the core consistency, and the number of iterations of the algorithm [19]. In this study, the core consistency method was employed. The core consistency represents the resemblance between the Tucker3 core and the PARAFAC core [14]. A Tucker3 model is similar to the PARAFAC model, and both of them are an extension of bilinear factor analysis to tensors [20]. The principle of the core consistency is as follows: (1) the PARAFAC model is available when the core consistency value is greater than 90%; (2) the PARAFAC model is not available when the core consistency value is less than 50%; (3) the PARAFAC model is probably available when the core consistency value is between 50 and 90%. The core consistency is defined as follows:=∑Core-Consistency1−^𝑅𝑡=1∑^𝑅_𝑗=1∑^𝑅_𝑘=1𝑔_𝑖𝑗𝑘−𝑡_𝑖𝑗𝑘²𝑅×100, (4) where g_𝑖𝑗𝑘 and 𝑡_𝑖𝑗𝑘 are the Tucker3 core and the PARAFAC core, respectively; 𝑅 is the factor number. In the PARAFAC core, 𝑡_𝑖𝑗𝑘=1 if 𝑖=𝑗=𝑘, otherwise 𝑡_𝑖𝑗𝑘=0, and in the Tucker3 core, 𝑔_𝑖𝑗𝑘 can be nonzero for all 𝑖, 𝑗, and 𝑘. This method for determining the factor number has been successfully applied to multiple neural data [13, 21].

2.4.4. Statistics

Data were analysed in the time, frequency, and channel domains for all subjects. For the channel domain, the average energy for each stimulus was calculated for all channels and compared between the two groups, with outlier detection based on the generalized extreme studentized deviate (GESD) [22]. The average energy distribution for each group was obtained by averaging across all subjects. Cortical locations were grouped into 5 zones (frontal (FPz, FP1, FP2, Fz, F3, F4, F7, F8), central (Fz, Cz, FC3, FC4, C3, C4, CP3, CP4), occipital (Pz, CP3, CP4, P3, P4, P7, P8, O1, O2), left temporal (F7, F3, Fz, FT7, FC3, T7, C3, TP7, CP3, P7), and right temporal (F8, F4, Fz, FT8, FC4, T8, C4, TP8, CP4, P8)), as shown in Figure 1(C). The active ratios (average energy of one zone/the sum of energy across all zones) were assessed using the 𝑡-test. In time and frequency domains, the differences between pain and control groups were evaluated by 𝑡-test.

The mean age of the chronic pain group (49±11 years, 𝑛=13) was significantly greater than the pain-free group (39±11 years, 𝑛=13) (𝑃<0.05, 𝑡-test). Multiple linear regression (MLR) was used to investigate the effect of age on the frequency domain parameter. The pain-free subjects and chronic pain subjects were represented by 0 and 1, respectively, and the effect of age and subject grouping differentiated by MLR as shown in Figure 4. In the MLR analysis, frequency was the dependent variable, while the age and subject grouping were independent variables. The effect of MLR was evaluated by the 𝐹-test.

3. Results

An example showing the epoch decomposition process using the PARAFAC method is illustrated in Figure 1(B). Wavelet decomposition was employed on the raw EEG recordings, generating time-frequency information corresponding to each channel. The frequency range for the wavelet decomposition was 1–50 Hz with an interval of 0.5 Hz. The factors were then extracted by the PARAFAC model from each epoch for every subject. Every factor was decomposed into three-way information at the time-frequency-channel domain. In the case shown (Figure 1(B) a, b, and c), the three-way information for one factor revealed a neural response frequency in the frontal region of 4 Hz and response time of 220 ms after the stimulus. Three-way information for each subject was obtained by averaging the results across all epochs. Comparisons between groups of the spatial topography, frequency, and time responses are shown in Figure 2.

Figure 2 

Results at the level of group corresponding to the four different stimuli. Left: the average energy distribution (pain, random pain, random touch, and touch, resp.). Middle: the statistical results in the frequency domain. Right: the statistical results in the time domain.

3.1. Comparison in Frequency Domain

As shown in Figure 2 (middle), the frequency response in the chronic pain group (around 4 Hz) was lower than that in control group (around 5.5 Hz) and was statistically significant for the pain stimulus (chronic pain group: 3.456±1.716Hz; control group: 5.608±2.315Hz; 𝑃<0.05, 𝑡-test).

3.2. Comparison in Time Domain

As shown in Figure 2 (right), the response time was around 250 ms in both groups and showed no significant difference for any of the stimuli.

3.3. Comparison in Location

As shown in Figure 2 (left), the topographical analysis indicated that the active zone was mainly in the frontal region in pain group and in the central region in the control group. These differences were quantified by comparing the active ratios in the 5 cortical regions (see Figure 1(C)). The results are shown in Figure 3(a); in the frontal region, the active ratio of neural response in the chronic pain group was significantly higher than that in the control group for all 4 tests stimuli (chronic pain group: 0.03551±0.00085, control group: 0.03459±0.00080, pain:𝑃<0.01; chronic pain group: 0.03557±0.00088, control group: 0.03469±0.00076, random pain:𝑃<0.01; chronic pain group: 0.03551±0.00077, control group: 0.03482±0.00074, random touch:𝑃<0.05; chronic pain group: 0.03545±0.00061, control group: 0.03485±0.00073, touch:𝑃<0.05, 𝑡-test). In comparison, in the central region, the active ratio of neural response in the chronic pain group was lower than that in the control group and was statistically significant for all but the touch stimulus (chronic pain group: 0.03555±0.00043, control group: 0.03624±0.00051, pain:𝑃<0.001; chronic pain group: 0.03557±0.00077, control group: 0.03608±0.00050, random pain:𝑃<0.05; chronic pain group: 0.0356±0.0006, control group: 0.03621±0.00043, random touch:𝑃<0.01; chronic pain group: 0.0355±0.00042, control group: 0.03591±0.00065, touch:𝑃=0.0526, 𝑡-test). In other regions, there were no significant differences between the two groups.

(a)

(b)

Figure 3 

(a) Statistical results of the active ratio between the two groups corresponding to 5 zones under the four different stimuli. (b) Within-group comparison of the active ratio between the frontal and the central zones in the chronic pain group (left) and the corresponding result in the control group (right).

Figure 4 

Influence of age and pain status on the frequency responses of subjects by MLR during the four different stimuli.

Furthermore, we were interested in the differences in active ratio between the frontal and central regions within each group. As shown in Figure 3(b), in the pain group, the active ratios were not significantly different, whereas in the control group, the active ratios in the central region were significantly higher than in the frontal region for all stimuli (frontal: 0.03459±0.00080, central: 0.03624±0.00051, pain:𝑃<0.0001; frontal: 0.03469±0.00076, central: 0.03608±0.00050, random pain:𝑃<0.0001; frontal: 0.03482±0.00074, central: 0.03621±0.00043, random touch:𝑃<0.0001; frontal: 0.03485±0.00073, central: 0.03591±0.00065, touch:𝑃<0.001, 𝑡-test).

3.4. Effects of the Age on the Time Frequency and Active Ratio

The mean age of the subjects in the pain group was significantly greater than in the control group. This study considered the effects of the age on the time, frequency, and active ratio parameters. Firstly, linear regression analysis was carried out on age versus response time, frequency, and active ratio; the only significant correlation was observed between age and response frequency. MLR was used to determine whether the difference in the frequency domain between two groups could be attributed to an age effect and is shown in Figure 4. The line marked by “∗” was the vertical mapping of the line of control group to the same plane with the line of the chronic pain group. From the comparisons of these lines, it is clear that the frequencies of the chronic pain group were significantly lower than those of the control group at the same age for the pain (𝑅²=0.4124, 𝐹=24.131, 𝑃=0.0017), random pain (𝑅²=0.3398, 𝐹=6.177, 𝑃=0.0069), random touch (𝑅²=0.4064, 𝐹=8.2145, 𝑃=0.0019), and touch (𝑅²=0.2156, 𝐹=3.2985, 𝑃=0.0542) stimuli. The differences in the frequency domain between the two groups therefore cannot be attributed entirely to the difference in age between groups.

4. Discussion

4.1. PARAFAC as an Effective Method for EEG Decomposition

In this study, the PARAFAC method was applied to somatosensory evoked potential recordings to analyze EEG time-frequency-channel domain [23] characteristics in chronic pain subjects. The advantage of this method is that it can extract more information in comparison with the two-way models (PCA and ICA) but also takes into account the frequency content of the signals in specific time periods across different channels [20]. The PARAFAC method has been successfully used to characterise the structure of epileptic seizure [14, 24, 25]. Analyzing the EEG signals in the time-frequency-channel domain, we found that the response latency was about 250 ms for all stimuli in both groups, while the frequency responses of the chronic pain group were lower than those of the control group. Furthermore, topographical analysis showed that the chronic pain group exhibited predominantly frontal cortical activity, compared to central activation in the control group. These findings are in accordance with previous results showing a response latency of about 250 ms [26], an increased size and reduction in population spike frequency [27], and restriction of some ERP components to frontal-central regions in patients with fibromyalgia syndrome [28]. Our findings indicate that the PARAFAC method is an effective technique for extracting the characteristics of multiple EEG recordings in the time-frequency-channel domain.

4.2. Activity Regions Involved in Chronic Pain

Previous studies have shown that the scale distribution of laser-evoked potentials (LEPs) around the chronic pain ERP components extends into vertex and frontocentral leads in fibromyalgia syndrome (FS) patients, indicating more widespread nociceptive activation outside the cortical hand area [28]. One question addressed in this study is whether the active ratio of frontal or central regions versus the whole cortex is different between the two groups. In this study, the active ratio was used to assess the activity within different regions of the cortex. The statistical results showed that the active ratios of the chronic pain group were significantly higher than those of the control group in the frontal region. The active ratios calculated from the other cortical regions showed no differences between the two groups. These results imply that evoked EEG activity in the frontal and central cortical regions may help discriminate between chronic pain and pain-free subjects. These findings are in agreement with the previous studies showing that some ERP components are more restricted to frontocentral regions in patients with fibromyalgia syndrome [28].

In keeping with the above findings, the active ratio in the central region was significantly higher than that in the frontal region in the control group, and vice versa for the chronic pain group. These results further indicate that the frontal cortical regions are involved in somatosensory processing in the chronic pain condition compared to central regions for pain-free subjects, in accordance with [29–32]. These findings suggest that passive functions (emotion, attention, etc.) are presented more frontally [33], and pain beliefs influence patients behavioral and psychological functioning because of their persistent pain experience [34].

4.3. Relationship between Frequency and Chronic Pain

In this study, we found that the cortical neural responses to somatosensory stimuli occurred around 250 ms in both of the groups and that the response frequency in the chronic pain group was lower than in the control group. In particular, the dominant activity was in the delta frequency range (around 4 Hz) for the chronic pain patients, compared to the theta frequency range (around 5.5 Hz) in the controls. Taken together, these results indicate that frontal cortical activation and a lower response frequency are characteristics of evoked EEG responses in chronic pain subjects.

In imaging studies, the functional connectivity between cortical structures receiving input arising from nociceptors has documented that experimental pain is processed in multiple pain-related areas, often characterized as a “pain network” [35, 36]. This “pain network” is not fixed but changes as a function of the pain-related task [37]. Moreover, in the study of cortical pathophysiology, pain-related neural networks are larger in patients with sympathetically mediated chronic pain (SMP) compared to acute pain states [32]. Therefore, the pain-related network of chronic pain patients is more extensive than that of the controls. This may be partly due to the effects of the patients’ persistent pain experience, resulting in processing experimental pain accompanied with subjective experience. The larger pain-related network is accompanied with a lower frequency rhythm, which is consistent with data from isolated hippocampal slice experiments showing a reduction in the frequency of population spikes with increasing size of the participating neuronal population [27].

Theta activity is associated with alertness, attention, and the efficient processing of cognitive and perceptual tasks [38], while delta band activity is associated with pathological conditions associated with impairment of brain networks such as amnesic mild cognitive impairment (MCI) and Alzheimer’s disease (AD) [39]. Chronic pain is also found to be frequently associated with psychiatric disorders [40]. Thus, the delta rhythm of chronic pain patients is suggestive of impairment to pain-related brain network processing.

4.4. Influence of Age on the Response Frequency

In this study, the mean age of the chronic pain group was significantly greater than that of the control group. MRL was applied to analyze the influences of age and pain status on the frequency responses of subjects. Our results showed that the frequency responses of the chronic pain group were significantly lower than those of the control group across all ages. Thus, the differences in the frequency domain between the two groups cannot be attributed entirely to the difference in age between groups.

In summary, the PARAFAC method with a continuous wavelet was used to extract time-frequency-channel domain information from somatosensory-evoked EEG recordings from chronic pain and pain-free subjects. We found that the response latency was about 250 ms, the chronic pain group had lower response frequency, and the central and frontal regions were the crucial regions of cortical activation. Further analysis indicated that the frontal regions were more involved in the chronic pain condition than the control condition. Application of MLR to the analysis of the relationship between frequency and age showed that the lower frequency response in the chronic pain group was not attributable to the difference in subject age. The conclusion from these findings is that the PARAFAC method is an effective tool for characterising multiple EEG recordings in the time-frequency-channel domain.

Acknowledgments

This research was partly supported by the National Natural Science Foundation of China (nos. 61025019, 90820016), the Program for New Century Excellent Talents in University (NECT-07-0735), and the Natural Science Foundation of Hebei, China.