Six subjects (4 males and 2 females, numbered 1–6) aged 22–25 years (standard deviation [SD] = 3.0) were recruited in this study. All participants were healthy, without any history of medical illness, neurological or psychopathological disorders, and none had a history of chronic pain. All the experimental protocols of this study were approved by the ethics committee of China Medical University. Moreover, all subjects were fully informed about the procedures, risk, and benefits of the study, and written informed consent was obtained from all subjects before the study.

The algorithm design and performance validation in this pain intensity measurement were implemented for the databases constructed from physiological signals of each subject. The signals were obtained with the procedure described in Section Physiological Signals and Preprocessing and Pain Induction Protocol. Moreover, we continuously collected 1-week data from the same participant, thus generating a 7-day dataset. For the database of each subject, we randomly chosen 75% of the samples as the training set, and the rest as the testing set. For the 7-day database, data from six randomly chosen days were used as the training set, and the rest were used as the testing set (leave-one-out method).

Acquisition of high-quality physiological data is of utmost importance for the pain intensity recognition system. The selection of available physiological signals that are used as input is the first consideration in the pain detection system. It is expected that the physiological signals reflect the effect of pain on the activity level of the nervous system. However, unlike self-reporting or facial recognition of pain, where the truth class labels of pain intensity for a given sample are self-evident, a high-confidence physiological signal in an underlying pain intensity is not easy to be acquired. Furthermore, it is difficult to determine whether the variability of the physiological signals is specific to pain or whether they represent a general response to factors such as thought, cognition, emotion and environment.

Since we were decided to exploit a practical quantitative and convenient system, there was a limitation on the selection of available signals. Though facial electromyograms and electroencephalograms would be expected to be useful, the attachment of electrodes to the face and scalp is complicated and not suitable for practical application. In our study, the selected physiological signals were blood volume pulse (BVP), electrocardiogram (ECG), and skin conductance level (SCL). These signals reflect the activity level of the autonomic nervous system, which is connected with the secretory activity of cardiac muscles and internal organs. As a wearable device, we use the Infiniti 3000A platform developed by Thought Technology Ltd., Quebec, Canada. The Infiniti 3000A apparatus is a compact sensing platform with 10 isolated channels for recording signals. It is small enough to attach to a portable computer and integrates to commercial sensors via custom cabling. The sampling frequency was fixed at 256 Hz for all the three channels.

The BVP signal is derived from a photoplethysmographic (PPG) sensor that monitors blood volume in capillaries and arteries by emitting an infrared light through the tissues. Vasomotor activity, which controls blood vessel diameter, is regulated by the sympathetic nervous system (Babchenko et al., 2001). Hence, changes in BVP amplitude reflect instantaneous sympathetic activation. Most PPG sensors can be placed anywhere on the body, with the finger as the most common location for recording a BVP signal. In this study, A BVP-Flex/Pro sensor was placed against the palmar surface of the middle finger of the right hand with an elastic strap or a short strip of adhesive tape to acquire BVP signals (Figure 1A). In the preprocessing phase, a 4-order Butterworth bandpass filter with gain 3, cutoff frequencies [30, 200] Hz was applied to eliminate the bursts in the BVP signal.

The ECG, an electro-physiological signal that is associated with the electrical activity of the sinuatrial node, reflects the cardiovascular activity. Additionally, ECG responses to external stimuli (such as pain stimuli and stress) can produce large variability in a given subject's physiological signal (Sriram et al., 2009). Therefore, we can employ ECG signal to extract universal information about pain state or intensity. The ECG signal was measured from both surface of upper arms with the two-electrode method based on lead I in our experiment (Figure 1A). For the preprocessing of the ECG trace collected by using non-invasive surface electrodes, a typical high-pass and moving-average filter (Milanesi et al., 2008; Sidek et al., 2014) were used to remove a low frequency baseline drift due to electrode contact noise, respiratory effects, and motion artifacts. Otherwise, we employed an adaptive linear interpolation algorithm to perform baseline correction before feature extraction.

The SCL is another signal that can easily be acquired from the body surface. It reveals variance in the electrical properties of the skin due to the secretory activity of the perspiratory glands. As it is directly regulated by the sympathetic nervous system, it is a good indicator of psychological or physiological arousal level due to external stimuli. Changes in the SCL signal induced by a painful stimulus were validated to be used as a tool to monitor nociceptive stimulation and pain (Harrison et al., 2006; Storm, 2008). The SCL signal was traced from two Ag/AgCl electrodes attached to the tip of the index and ring fingers on the palm-side of the right hand in this study (Figure 1A). In the preprocessing phase, a moving average smoothing was carried out to subtract artifacts. Furthermore, the raw SCL data was down-sampled by a factor of 1/2 to reduce the computational complexity in feature extraction phase.

As in many pain researches, there are variety ways that can be used to induce pain, such as thermal or cold pain stimulation (Appelhans and Luecken, 2008; Kachele et al., 2016), mechanical pain elicitation (Matsunaga et al., 2005; Shankar et al., 2009), and electrical shock (Oliveira et al., 2012; Zhang et al., 2016). In our study, pain was induced by an electrical stimulator (MotionStim8; Medel GmbH, Hamburg, Germany), which can generate a current square wave with a certain pulse width. The amplitude and frequency of the current are adjustable. The value of the pulse amplitude can be used as an objective index of pain intensity. In the experiment, the frequency of MotionStim8 was set to 2 Hz. In order to avoid interference of the stimulator with the sensors, especially the ECG sensor, the stimulation electrodes were placed on the tibialis anterior muscle of the right leg, as far away as possible from the sensors (Figure 1A).

Before stimulation phase, the current intensity was calibrated according to the subjects' self-reports. An electrical stimulus with 20 mA induced the sensation of pain was as pain starts (threshold). An electrical stimulus with 40 mA induced intense pain was as barely bearable pain (tolerance). Then, we divided the range between threshold and tolerance into 3 equally spaced intervals, 20-mA stimulus (stim20, L1), 30-mA stimulus (stim30, L2), and 40-mA stimulus (stim40, L3) with additional pre-stimulus (baseline, L0). During the stimulation experiment, each of the different current intensity was proceed for 1 min followed by a recovery period of 1.5–2 min (Figure 1B). To eliminate stimulation adaptability and time correlations, the sequence of stimuli and duration of the recovery period were randomly designed. Pain elicitation was executed 30 sessions for each of the 3 calibrated intensities (L1–L3) in the same day. With the addition of baseline phase, a total number of 120 stimulation trials were obtained for each subject. Each session took approximately 15 min. To minimize motion artifacts, the participants were requested to be as relaxed as possible during the stimulation phase. Moreover, the recording scenario that led to the multi-day dataset was continuously carried out 7 days for one person in the same condition. Figure 2 shows physiological traces from a participant in the four pain states.

From the Figure 2, we can see that the amplitude of physiological signals obviously varied according to the level of electrical stimulus. Based on this observation, we gathered a variety of basic statistical features to evaluate the response to the pain-inducing stimuli, such as maximum, minimum, median values of data. Furthermore, these features can easily be calculated in an online way, which makes them more suitable for real-time recognition system of pain intensity.

The statistical features can be calculated for each of the signal channels as follows. Let the three preprocessed signals—the digitized BVP, ECG and SCL waveforms from any one of the four pain states segments—be designated by Xⁱ, i = 1, 2, 3 respectively. Let Xni represent the value of the n^th sample of the i^th signal, where n = 1, 2, …, N with N in range of 3,840 to 7,680. Let X~ni refer to the normalized signal (zero mean, unit variance), formulated as:

where μⁱ and σⁱ are the means and standard deviations of Xⁱ as presented below. Let fji represent the j^th features of the i^th signal.

• The means of the signals: f1i:

• The standard deviations of the signals: f2i:

• The means of the absolute values of the first differences of the signals: f3i:

• The means of the absolute values of the first differences of the normalized signals: f4i:

• The means of the absolute values of the second differences of the signals: f5i:

• The means of the absolute values of the second differences of the normalized signals: f6i:

• The minimum: f7i:

• The maximum: f8i:

• The minimum ratio: f9i:

Here, N represents the length of the signal.

• The maximum ratio: f10i:

Here, N represents the length of the signal.

• The range of the signal: f11i:

Besides, we also calculated the median value of the signal as f12i. Hence, each channel was characterized by 12 features. Considering the time lag of pain reaction, we selected 30 s signal for feature analysis in the 1-min data. For obtaining the subtle changes of the signals in a trial, we segmented the 30 s signal by a sliding window with 3 s sliding step (the window length was 3 s) to obtain the feature windows. Hence, a sample dataset of 40 × 36 can be obtained for each experiment session consisted of 4 trials. For single subject with 30 sessions, the size of sample dataset was 1,200 with 36 feature dimensions.

The BVP, ECG, and SCL are dependent on each participant's initial physiological level. Even when these signals are measured from the same person, they are likely day-dependent due to variations in mental state affected by emotion, variations in physiology caused by sleep or diet, or variations in the sensor's connectivity with skin (Sun et al., 2010). To eliminate the intra-subject and day-to-day variation factor, a normalization stage for each feature was applied in feature vector set F. Refer to (13), the first was to subtract the minimum value from each feature. Then, the features are divided by the range to make all the values lie between 0 and 1.

where, V_min(m) = min{F_n(m)}, V_max(m) = max{F_n(m)}, ∀n ∈ |F|.

Note that not all of these features are independent. Some features are a linear combination of other features (e.g., θ~1i and θ~2i) or a non-linear property. To obtain high-quality features with significant discriminatory information, a hybrid combination of genetic algorithm (GA) and principal component analysis (PCA) was proposed. This combination was motivated by the roles of the two algorithms: GA selects from a feature set to remove useless and redundant features, while PCA linearly transforms a feature set. Since feature selection is non-linear, the hybrid of these two ways could provide a more powerful separation than either alone.

The GA is a biologically inspired stochastic algorithm based on some optimization criterion that mimics natural evolution. The GA is naturally applicable to feature selection since the problem has an exponential search space (Oh et al., 2004). Furthermore, the distinct specialty of this algorithm is that it always maintains a set of solutions (called chromosomes or individuals) in a population. To simulate the biological evolution, the GA selects fitter chromosomes at each generation by some genetic operations, such as crossover and mutation. The traditional procedures of steady-state GA are described below.

For the feature selection problem, a finite bit string with d binary digits is usually used to represent chromosomes. A binary digit represents a feature. A value of 1 or 0 in the chromosomes means that the corresponding feature is selected or removed, respectively. For example, a chromosome 10100010 means that the first, third and seventh features are selected. Since we wanted to improve recognition rate, the classification accuracy from the classifiers was directly used as the fitness function to evaluate a chromosome (Lin et al., 2015). Based on the initial fitness evaluation, the chromosome for the next generation was selected by a selection mechanism that ensures fitter chromosomes have a higher probability to survive. In our design, a traditional roulette wheel selection strategy was applied. We also used the standard two-point crossover operator to obtain the offspring and Gaussian mutation operations at probability on each of the offspring produced from crossover. The crossover rate and mutation rate were set to 0.8 and 0.01, respectively. Besides, the size of population and the maximum generation in GA were set to 100 and 30, respectively. These parameters are tunable to be suitable for a specific dataset to improve performance.

The principal component analysis (PCA) is a well-known statistics method of multi-variable analysis and used orthogonal transformation to convert a set of observations of possibly correlated variables into a set of values of linearly uncorrelated variables (Karamizadeh et al., 2013), also known as principal components. Moreover, PCA is a tool to reduce dimensionality by searching a linear projection matrix of the dataset to a feature space with lower dimensions while retaining most of the information. This linear projection matrix is a meaningful basis to filter out the noise and reveal hidden information in the original feature space. The key advantages of the PCA are the decreased requirements for memory and capacity, it low noise sensitivity, and increased efficiency in a lower dimension. By solving the covariance matrix, a series of eigenvalues with the corresponding eigenvectors are computed and rank-ordered descendingly by quantifying how “principal” each variance is. The k largest eigenvalues are selected and the corresponding eigenvectors (k principal components) are used to construct the linear projection matrix. Then, a feature space with N dimensions is transformed to a space of k dimensions by this reduction matrix. For a detailed description of the PCA algorithm about mathematical derivation, the reader is referred to Bro and Smilde (2014).

A hybrid GA with PCA method is proposed, evaluated and implemented for comparison. As mentioned above, the GA algorithm provides a subset of features selected from raw features. Instead of feeding the PCA algorithm with all possible features, we used the selected subset of features as the input to the PCA. Note that the PCA algorithm is applied here for transforming features into a linearly uncorrelated space.

Previous studies have performed qualitative analyses of pain based on a single physiological signal. For example, Harrison et al. revealed SCL as a promising indicator to predict the presence or absence of pain in hospitalized infants during painful and non-painful medical procedures (Harrison et al., 2006). However, in this study, we applied pattern recognition techniques to quantify pain intensity based on a set of physiological signals. Since there was not yet a priori knowledge of pain recognition, we originally established a quantitative model for pain intensity from individual's single physiological signal. For each signal, we extracted 12 statistical features to construct three types of feature set, including BVP-features, ECG-features, and SCL-features. The classification accuracy of the LDA algorithm was evaluated by using these three types of feature set. Figure 4 presents the recognition results gained with LDA algorithm for each physiological signal in each of the six study subjects.

From the Figure 4, the experimental results provide evidence to support our hypothesis that a single signal can be used to quantify pain intensity in some degree, with the average classification accuracy of 50% above. Furthermore, since the BVP and SCL signal were highly affected by the intensity of pain stimuli in our experiment, the classification accuracies were even higher than that of ECG signal alone. Especially, the average classification accuracy of SCL signal was 68.39% (±2.32%), and that of ECG signal was 53.01% (±3.40%). Unlike ECG, it can be observed that even the original SCL signal already presents an obvious variance with the pain stimulation in Figure 2. As the experimental results suggest, the extracted SCL features are a good indicator to identify the pain intensity. This is particularly suitable when it comes to integrate the proposed methods into wearable device. Since the SCL sensors are less intrusive and more convenient attached than BVP or ECG sensors, it can be easily embedded into wearable devices. Although the performance of LDA model for a single signal is considerable, the correct classification rate is relatively low. For improving the classification performance and obtaining better models, we considered applying multiple physiological signals with all-feature combination to quantify pain intensity. The idea was inspired to take full advantages of pain information contained in various signals.

Since we used multiple signals (BVP, ECG, and SCL), the total number of extracted features per sample was 36 (12 features per signal). We also selected the LDA algorithm as the classifier. The recognition accuracies of running 10 trials LDA are shown in Figure 5. The average accuracy of the LDA model was 42.49% (±14.31%) for subject 1, 53.78% (±20.84%) for subject 2, 65.12% (±16.57%) for subject 3, 51.07% (±20.52%) for subject 4, 54.30% (±13.62%) for subject 5, and 55.52% (±11.43%) for subject 6, respectively. The maximum classification accuracy of 10 trials LDA algorithm in the condition of multiple signals was approximately 90% (such as subjects 2 and 3). However, the variation of classification accuracy was large, ranging from 19.64% (subject 1) to 87.61% (subject 3). This finding indicated that the robustness of the multi-signal model established by LDA was insufficient. This may be because some extracted features are not independent and may be correlated. Furthermore, some features do not contribute or contribute little to the performance of the classification. A poor feature combination may greatly degrade the performance of the classifier. Hence, feature selection was needed to select a subset of good features.

In this study, we applied the genetic algorithm technique to seek and identify the potential informative features combinations, which contribute most to the performance of the classification. From other researches, it is obvious that the best accuracy of the LDA by using cross-validation method is considered as the indicator of the fitness function. Table 1 provide the experimental results of GA selection. The GA was run 15 times; each row lists the results from the corresponding run. In Table 1, the column “Best generation” indicates the generation number in which the optimal set of features was selected. And the average number of selected features was roughly 18. The performance of LDA was significantly improved with average accuracy of 72% above. In detail, the feature f2i, f3i, f4i, f5i, f6i, and f7i computed from the three physiological signals were almost repeatedly selected by each running of the GA, yielding 18 features: μⁱ, σⁱ, θ1i, θ~1i, θ2i, θ~2i, i ∈ (BVP, ECG, SCL). It means that these features are the more robust feature combination than those that are rarely selected in the context of this pain experiment.

Although the trend of the feature combination proves their efficacy in distinguishing pain intensity, the 18 features selected by GA are not the best informative features for classification. This is because the genetic algorithm cannot easily eliminate the linear features. To solve this problem, we applied the feature reduction technique. Figure 6 shows samples of 18-dimension feature space projected onto the first three features from subject 1. The selected feature vectors from the same class aggregated a cluster with a large amount of variation, whereas the feature vectors from the different classes overlapped dramatically. Then, a PCA method was applied to reduce the redundant linear features. After data transformation, the feature vector samples within same class were clustered more compactly than did samples before transformation, and the feature vector samples from different classes are farther away from one another (see in Figure 7).

After the feature processing stage by using the hybrid GA with PCA, the samples were fed to the LDA classifier to investigate the performance. In addition, we also tested the KNN and SVM algorithms for comparison. The corresponding recognition accuracy is presented in Figure 8. Using the optimized feature set, the average classification accuracy of three classifiers reached 80% or higher. Particularly, the performance of LDA algorithm was better than the other two algorithms, with 98% average accuracy. Those results prove that the proposed feature processing method is effective. And, the LDA classifier is more suitable in physiology-based pain intensity recognition for a single-subject scenario.

In addition to discussing the recognition model of a single subject, we also analyzed samples constructed from multiple subjects. A sample set of multiple subjects was generated based on the six participants. Similarly, GA-based feature selection and PCA-based dimension reduction technique were carried out on the sample set comprising multiple subjects. The LDA, KNN, and SVM algorithms were implemented to train the corresponding models. Figure 9 shows the pain intensity recognition accuracy of running 10 times for the multi-subject model. The average categorize accuracy was 84.28, 83.94, and 96.47% for LDA, KNN and SVM, respectively. Overall, these three algorithms effectively identified pain intensity. The relevant multi-subject models established by LDA, KNN and SVM are feasible and appropriate. Besides, the multi-subjects model gained from SVM algorithm was better than the other two models. Every time we used the training SVM model to test the rest samples, the predicted accuracy was 92% or higher. Therefore, the trained SVM model can be used as a general model to quantify pain intensity for a multiple-subject scenario.

Meanwhile, between-subjects' experimental validation was conducted by using the LDA, KNN, and SVM in a leave-one-subject-out cross-validation method, which is standard when dealing with multiple subjects. Here, the recognition modes were trained with dataset from randomly selected six subjects to predict the pain intensity of the remaining subject. It was repeated until all of the subjects had been part of the testing dataset. For a given between-subject scenario, we repeated 20 times and calculated the average accuracy for each pain intensity, respectively (see Figure 10). In detail, the average recognition accuracy of pain intensity with baseline, stim20, stim30, stim40 was 78.01, 90.10, 82.35, 82.90%, respectively. This is a more realistic condition in which the distribution of the training dataset and the testing dataset are dissimilar due to inter-individual variability. Moreover, the SVM algorithm was performed better than the KNN and LDA, with average accuracy of 91.18, 76.14, and 82.83%, respectively. These results demonstrate that the SVM classifier outperforms the other two classifiers for a between-subject scenario. Furthermore, the proposed system in the context of pain prediction would mitigate the inter-personal variance while preserving precise scalability.

Furthermore, we also conducted experimental evaluation by using the three types of classifiers in a multi-day scenario. For one person, we continuously carried out the pain induction experiment of 7 days and collected a seven-day dataset. A leave-one-day-out cross-validation method was applied to obtain the pain intensity of one person. The results of multi-day model are given in Figure 11. The average recognition rate was 71.02, 81.39, and 76.93% for KNN, SVM, and LDA, respectively. It implied that the trained SVM model seemed to be more suitable as a special model to quantify pain intensity of multiple days for one subject. In addition, the overall accuracy was roughly reduced from 87% of 1-day scenario to 70% of multi-day scenario. This may be mainly because that the signals measured from the same person are likely day-dependent. The quality of signals is affected by the intro-individual variance in mental state and environmental factor. However, our proposed pain intensity system is relatively fitted for practical applications of continuous and long-term pain monitoring.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.