2.1.

Subjects

Eighteen right-handed volunteers (three females) were considered in the study. All participants were of Chinese ethnicity and dark brown hair color, $\text{mean age}\pm \text{standard deviation}$ ($31.9\pm 5.5$). No participants reported a prior history of neurological or psychiatric disorder, a current unstable medical condition, or under medication at the time of testing. All participants were right-handed to avoid any variation in functional response due to lateralization of brain function. The experimental procedure was explained to the participants, and they had the opportunity to stop the procedure at any time if it was needed. Written informed consent was obtained from all volunteers after providing an explanation of the protocols in the study. Procedures and methods for this study followed the guidelines accepted by the Declaration of Helsinki. This research study was approved by the Taipei Medical University (TMU) and full-board review process of the TMU-Joint Institutional Review Board under contract number 201307010.

2.2.

Experimental Paradigms

Thermal pain perceptions (threshold and tolerance) in healthy volunteers were investigated following the standard procedures of the QST protocol.²² We defined pain threshold (low pain) as the least stimulus intensity at which stimulation becomes painful and pain tolerance (high pain) as the highest intensity of pain at which stimulus becomes unbearable. Both thermal threshold and thermal tolerance were obtained by heat and cold stimulation, and pain measurements were obtained on the dorsum of the left hand (nondominant hand). All thermal tests were performed using a pathway CHEPS (Medoc, Israel), this sensory analyzer delivers heat and cold to the skin with a thermode. Attached to the central unit is a hand-held controller that allows subjects to press a button to stop the heating/cooling of the thermode. The subjects were seated on a comfortable armchair with the left arm resting on a table. Previous studies have reported no significant difference in QST between the right and left sides of the body.²²

The pain experiment involved applying cold and heat to the skin to induce pain. In this test, subjects are exposed to gradually increasing or decreasing temperatures with a thermode and they pressed a button when they experienced pain (threshold test) and highest intensity of pain (tolerance test). The temperature of the thermode, just as it becomes painful or unbearable, is the thermal pain threshold or thermal pain tolerance, respectively. The pain threshold and tolerance for each subject were estimated as the average from three consecutive trials of cold and heat stimulation. Figure 1 summarizes the stimulation paradigm and is explained below.

Fig. 1

Stimulation paradigm. Pain threshold test was first measured and followed by pain tolerance test. In each test, cold and heat stimulus were applied on the back of the hand of each subject.

Similar to Rolke et al.,²² our protocol was divided in two tests: the thermal pain threshold (low pain) and the thermal pain tolerance (high pain), with a 2-min rest between both tests. After an initial 60-s rest, the thermal pain threshold was first measured by three consecutive trials of cold stimulation, followed by 60-s rest and then, three consecutive measurements of heat stimulation. Thermal pain tolerance was measured by following the same procedure for thermal pain threshold. In both tests, there was a 30-s interval between each trial. The three consecutive measurements from each test were averaged, and the block-averaged signal ($\sim 30\mathrm{s}$) from each thermal test was used for the classification task. Thermal measurements were obtained with ramped stimuli of $1°\mathrm{C}/\mathrm{s}$ that were terminated when the subjects pressed a button, then, the temperature returned to baseline (32°C) for the next stimulation with a rate of $10°\mathrm{C}/\mathrm{s}$. The probe has a safety cut-off to prevent any harm to the skin. The low cut-off temperature was 0°C and the high cut-off temperature was 50°C. The contact area of the thermode was $9{\mathrm{cm}}^2$. The testing took place in a quiet and temperature-controlled (20°C to 22°C) room at TMU. All experiments were carried out by the same researcher to reduce variabilities between samples. Based on the threshold and tolerance of cold and heat stimuli, the fNIRS data were organized into four categories: (1) low-cold (low pain), (2) low-heat (low pain), (3) high-cold (high pain), and (4) high-heat (high pain). These categories (1 to 4) were used to label the database and employed as classes for the classification task.

2.3.

fNIRS Recording, Preprocessing, and Channel Selection

fNIRS recordings were acquired using a multichannel optical topography system, Hitachi ETG-4000 (Hitachi Medical Corporation, Japan). This equipment uses two wavelengths of near-infrared light to measure changes of hemoglobin concentration in the cerebral cortex; HbO at 695 nm and HbR at 830 nm.⁹ The spectrometer is equipped with a 24-channel cap, configured in 12 channels (rectangles) per hemisphere (Fig. 2). Each hemispheric probe consists of five sources (red circles) and four detectors (blue circles) that provide 12 source–detector pairs. The head probe has a predefined source–detector separation (SDS) of 3 cm (left panel Fig. 2) and an NIR-light head penetration of about half the SDS distance.²³ According to the EEG 10 to 20 system, the measuring probes were centered on the C3 and C4 positions. In this study, we used only the HbO signals because HbO signals exhibit a better signal-to-noise ratio than HbR signals.^24,25 The sampling rate was 10 Hz.

Fig. 2

Channel location and configuration. Channel probes were located around the C3 and C4 areas. For classification purposes, two sets of channels were tested, all 24 channels and 8 channels from a RoI. The selected channels from the RoI are four channels on the right hemisphere (Ch4, Ch6, Ch7, Ch9) and four channels on the left hemisphere (Ch16, Ch18, Ch19, Ch21), which are positioned over the primary somatosensory region. Source–detector distance was 3 cm.

Raw fNIRS data are normally contaminated by different sources of noise, and preprocessing was required. To minimize effects of physiological components (e.g., cardiac, respiration) and electrical noise in fNIRS data, a finite-impulse forth-order low-pass filter with a cut-off frequency of 0.1 Hz was used. Then, motion artifacts were removed by following a denoising procedure similar to Molavi et al.,²⁶ which uses discrete wavelet coefficients to find and eliminate those coefficients (outliers) that do not belong to the probability distribution and therefore behave as motion artifacts; four subjects presented motion artifacts. Raw signals were converted to changes in HbO and HbR concentrations using the modified Beer–Lambert law applying literature differential path length factor values of 6.51 for the 695 nm and 5.86 for the 830 nm light sources.²⁷ In order to reduce extracerebral hemodynamics (scalp) and systemic variables²³ (e.g., blood pressure, breathing rate, autonomic nervous system activity) existed in our data, we followed the method presented by Zhang et al.²⁸ we applied principal component analysis (PCA) to the baseline data, assuming that baseline data primarily contain the spatial patterns of systemic interference and identified the first PC (which accounted for approximately 90% of the variance); then, we attenuated the interference in the stimulus data by projecting the stimulus data onto the orthogonal subspace of identified spatial eigenvectors to obtain “clean” stimulus data. Finally, we normalized all fNIRS data to zero mean and unit standard deviation to avoid different baselines and amplitude variations.

The number of fNIRS channels used in this study was evaluated during the classification task and two sets of channels were selected for that purpose. The first set of channels is formed of all the 24 channels available, whereas the second set is focused on a specific region-of-interest (RoI) composed of eight channels. The reasoning behind this evaluation is that the fNIRS head probe was covering a wide area of the cerebral cortex, including the frontal and parietal lobes on both hemispheres. However, for this particular study, we were interested in evoked signals related to pain in the left hand, therefore, our RoI was the somatosensory region in the parietal lobe near the central sulcus on both hemispheres. Based on this hypothesis, we decided to use the channels around this cortical area and compare the performance of RoI-based channels and all the 24 channels. The RoI-based channels from both hemispheres were: Ch4, Ch6, Ch7, Ch9, Ch16, Ch18, Ch19, Ch21.

2.4.

Bag-of-Words Representation

This study is based on the BoW approach for document classification adapted for time series. BoW representations are very popular in text mining, document classification, texture analysis, scene classification, and robot localization due to its simplicity and good performance.^29,30 BoW has become well-established in the field of computer vision, where it is also referred to as bag-of-features. Recently, it has been adapted for the classification of time series in audio and speech recognition,^31–33 electroencephalogram (EEG) and electrocardiogram (ECG) signals,³⁴ and time series similarity.³⁵ In these works, time series are treated as text documents and sections extracted from the time series as words.

In the context of fNIRS time series representation, BoW can adopt a similar model to the approach used for image classification. It can be generalized in three steps: first, the detection of keypoints (features) in the fNIRS time series using sliding windows; second, keypoints are represented into words (codebook generation) using a quantization technique (e.g., $k$-means); and third, fNIRS time series are characterized by histograms of number of occurrences for each word observed in each time series. This histogram representation is used as input to a classifier (e.g., SVM) that finds boundaries between classes. These steps are exhibited in Fig. 3 and are described below.

Fig. 3

fNIRS time series representation: (I) sliding a window over the entire time series, (II) local segments are converted to feature vectors, local descriptors (e.g., wavelet) can be used to represent feature vectors, (III) vector quantization using $k$-means clustering, (IV) codebook (or dictionary) generation, each cluster center produced by $k$-means becomes a word (${\mathrm{W}}_k$), the codebook is universal for training and testing data, (V) each feature vector is assigned to nearest word in the dictionary, (VI) histogram representation containing total number of appearances of each word in the time series. For new time series (testing data), repeat steps (I), (II), (V), and (VI).

2.5.

Classifiers

The classification task is to assign a vector $X_i$ of input $\mathbf{X}=[X_1,X_2,\dots ,X_i]$ to one of the $j$ classes $Y_j$. For that purpose, a classifier has to learn a model from training data ${\mathbf{X}}_{\text{train}}$ (a dataset reserved for this purpose). We then use this model to make class assignments on a different dataset ${\mathbf{X}}_{\text{test}}$ (independent from training data). This model captures the relationship between features and class labels in the training data. The end goal of a classification problem is to make generalized predictions on unlabeled data from the same distribution using the trained model. Properly, a classifier can predict the label $Y$ of a given feature vector $X$ as a function $f$ of the type $Y=f(X)$.³⁹ To obtain the best model on our pain database, two state-of-the-art algorithms were compared, the K-nearest neighbor (K-NN) and SVMs for multiclass classification (one-versus-one voting).

The K-nearest neighbor (K-NN) is one of the simplest classification algorithms for supervised learning and a robust classifier for time-series classification.⁴⁰ It uses the K features in the training set ${\mathbf{X}}_{\text{train}}$ closest (nearest) in feature space to unlabeled data $x_0$ to predict its output $y_0$. The algorithm basically uses a majority vote between the $K$ nearest training points and the query point. Similarity is defined based on a distance metric between two data points, for that purpose, we used Euclidean distance due to its simplicity and low computational cost.⁴¹ Different number of nearest neighbors ($K=\mathrm{1,3},\mathrm{5,7},9$) was also tested.

SVMs have been previously used in classification problems in the neuroscience field.^18,21 Its basic idea is to find a linear separating boundary with the maximal margin, and new data points are then classified according to which side of the decision boundary they fall.⁴² However, data in real-world scenarios are not easy to separate by a linear decision boundary; in these cases, SVM solves this problem by kernel functions to map the observations into higher-dimensional space in which the data can be separable. We used the LibSVM library⁴³ parametrized for a multiclass classification and the linear kernel was used for this study. One of the reasons to use a linear kernel compared to a higher-order kernel (e.g., radial base kernel) is that a linear kernel is faster, has less parameters to optimize, and it has been used in previous studies on pain research using brain imaging techniques.¹⁸

2.6.

Evaluation of Classification Models

The pain database was randomly split (on the subject level) in two parts, training and testing datasets. The training set was formed with 13 subjects and the testing set with the remaining 5 subjects. Therefore, the data are divided into two independent sets, one subset is used to train the classifier and another subset is used to test the performance of the classifier. During the training process, leave-one-out cross validation on the subject level is computed. One data point (1 subject) is held as the validation set while the remaining data (12 subjects) are used for training; the process is repeated 12 times, validating on a different subject in each iteration. During this process, the BoW parameters (step size, window length, and codebook size) were evaluated with a 1-NN and a linear SVM. The trained models were then tested using the remaining five subjects and the estimates (accuracy, sensitivity, and specificity) of predictive performance are then averaged across all the testing groups and reported as the final result.

3.1.

Thermal Threshold and Tolerance

Thermal threshold and tolerance of pain perception were obtained following a similar approach to the thermal test in the QST. By determining the QST thermal tests, we aimed to minimize the subjective nature of self-reported pain scores. In this way, the obtained measurements are based on temperature readings to classify a measurement as low pain (threshold) or high pain (tolerance) of cold and heat stimulation and not on self-reports. The measured (averaged) values obtained from each experiment are shown in Fig. 5.

Fig. 5

Thermal threshold and tolerance levels experienced by the participants (a) after cold and (b) heat stimuli. Horizontal red lines are the median values across all participants for each test.

The two plots show the threshold and tolerance temperatures of cold (Fig. 5, left panel) and heat (Fig. 5, right panel) stimuli across all participants. The first three measurements in each plot refer to the pain threshold while the last three measurements refer to the pain tolerance. The median ($\pm \mathrm{SD}$) temperatures (horizontal red lines in Fig. 5) in which participants first perceived pain (pain threshold) from cold ($12.45\pm 1.97,12.05\pm 1.93,12.45\pm {2.22}^{°}\mathrm{C}$) and heat ($42.70\pm 2.44,42.80\pm 2.75,42.95\pm {2.92}^{°}\mathrm{C}$) clearly presented differences from the highest intensity of pain (pain tolerance) the participants could take from cold ($3.40\pm 2.07,2.40\pm 2.71,1.45\pm {2.12}^{°}\mathrm{C}$) and heat ($48.00\pm 1.92,49.10\pm 2.11,49.70\pm {1.82}^{°}\mathrm{C}$). Median values (horizontal red lines) of cold pain threshold exhibited a statistically significant ($p<0.05$) higher response (Wilcoxon signed rank test, $p=0.0018$) than the median response to cold pain tolerance; similarly, median values of hot pain threshold showed a significant higher response (Wilcoxon signed rank test, $p=0.002$) to the average response to hot pain tolerance.

Habituation of noxious thermal stimulation was observed in all four experiments. As we can see in Fig. 5, group temperatures of pain threshold and pain tolerance slightly decreased and increased after each cold and heat stimulation, respectively. Therefore, we used Friedman’s test and looked for statistically significant differences between the three consecutive measurements (test 1, test 2, test 3) from each thermal test (e.g., cold threshold). However, the Friedman’s test did not show any statistically significant differences ($p<0.05$) between repetitions within cold threshold ($p=0.18$), heat threshold ($p=0.57$), cold tolerance ($p=0.22$), and heat tolerance ($p=0.13$) tests.

3.2.

Bag-of-Words Parameters

As a validity check on the effectiveness of the chosen parameters, the main parameters to build the BoW implementation were examined. Different step sizes, length of local segments, and codebook sizes were evaluated during the training process. The parameter evaluation included two classifiers (1-NN and linear SVM) and two channel sets (24 channels and 8 channels). The results of parameter tuning of the BoW representation are presented in Fig. 6 and are described below.

Fig. 6

BoW parameter tuning. (a) Segment (window) length, (b) codebook size, and (c) step size were evaluated; their corresponding classification performance while testing the parameters is also presented. Note that the units of segment length and step size are samples, and the units of codebook size are words.

We first evaluated the length of the sliding window by varying its size among the range of 32 to 352 in increments of 32 samples. In this process, we fixed the codebook size to 1000 and the step size to 4; these reference values have shown good results in similar studies.^30,34 The best results were obtained using a sliding window between 192 and 352 samples for the 1-NN classifiers and between 288 and 352 samples for the SVM classifiers. In all cases, the accuracy of the classifiers drops significantly with a length smaller than 160 samples. Following these results, lengths of 224 and 192 samples were selected for the 1-NN classifiers with 24 channels and 8 channels, respectively; lengths of 288 and 320 samples were chosen for the SVM classifiers with 24 channels and 8 channels, respectively.

Then, the influence of the codebook size on the classification performance was evaluated. We fixed the step size to 4 and used the best window length obtained in the previous step for each classifier. The results in performance during the training process are relatively stable with a codebook size between 50 and 1500 words for the 1-NN classifiers and between 500 and 300 words for both SVMs. It is important to reduce computational time and memory usage while having good discrimination, thus the best values for codebook size were 500 words for both 1-NNs, and 2000 and 2500 words for the SVM classifiers with 8 channels and 24 channels, respectively.

Finally, the last parameter to evaluate was the step size of the sliding windows. In this case, the best values for the windows length and codebook size were used to test different step size from one to nine samples. The evaluation showed stable results for most values, however, the best results were obtained with step size of four and six samples for the 24-channel and 8-channels 1-NN, respectively, and five samples for both SVMs. Therefore, the best BoW parameters were used to build the codebook, which was universal for both training and testing of each classifier.

In addition, an example of cluster centers obtained by $k$-means is presented in Fig. 7. The image presents the DWT’s approximation coefficients derived from the sliding window process. In this set of samples, we can observe distinctive patterns, i.e., increasing, decreasing, decreasing, and then increasing. Each of these centroids symbolizes the most representative value from a group of fNIRS segments with distinctive pattern similarities.

Fig. 7

Samples of cluster centers constructed by $k$-means on the zero-mean and unit-variance fNIRS using the DWT coefficients.

3.3.

Classification of fNIRS Signals in Response to Thermal Pain

The main goal of the classification task is to design a classifier that can make predictions using new observations. As mentioned before, we have four types of observations (classes); these observations refer to temperature level (hot and cold) and pain intensity (low and high). Based on these observations we trained our classifier to predict new observations to be either low-cold pain, low-heat pain, high-cold pain, or high-heat pain. With this in mind, we used the K-NN (with $K=1$, 1-NN) and the linear SVM algorithms with 24 channels and 8 channels (fNIRS channels) to design four classifiers able to estimate the type of pain and its intensity by using the BoW representation. Using the best BoW parameters for each classifier, we tested each classification model with unseen subjects ($n=5$), the results were averaged and summarized in Table 1.

Table 1

Performance of the K-NN (K=1) and SVM (linear) classifiers using the best results from parameter tuning of the BoW representation. The evaluation also includes two sets of channels, all 24 channels (Ch) and 8 channels, from the RoI. “Low” refers to pain threshold, while “high” refers to pain tolerance.

The overall performance of each classifier showed good generalization results with both channel sets. The best results were obtained by the 1-NN classifier with 92.08% and 91.53% using the 24 channels and 8 channels, respectively. On the other hand, the SVM classification models showed a slightly lower performance compared to the 1-NN classifiers, 91.25% using 24 channels, and 90.83% using the 8-channel set. In addition, we also evaluated different NNs with $K=3$, 5, 7, 9; however, in both channel sets, the 1-NN showed the best results.

4.1.

Building the Classifiers

The BoW representation showed its potential to characterize fNIRS signals, however, parameter tuning of such representation is important to obtain a trade-off between discrimination and generalization. While the step size of the sliding window showed stable results with different numbers of steps, the segment length and the codebook size presented more variability. Local information was extracted by continuously sliding a window along an fNIRS signal and BoW allows incorporating this information that appears at different temporal locations and with different scales in an efficient way. If the window segment is short, the local information loses meaningfulness due to the lack of data points for feature extraction and a large number of local segments will be extracted from each fNIRS signal, which increases the computational complexity due to the large number of segments to cluster. However, if the window segment is large, only a few local segments are obtained, which restricts the discriminative ability of the BoW representation due to the limited number of samples. On the other hand, the codebook size showed good results with medium-size codebooks, however, it varied for K-NN and SVM classifiers. In previous works on classification in biomedical signals, such as EEG and ECG^30,34,36 using BoW, the optimal codebook size has ranged from 25 to 8024 words, and this diversity is mainly due to differences between the nature and data type of each study. For our particular study, a compact codebook ($<500$ words) produced limited discriminative results because segments with a different profile can be assigned to the same cluster, whereas a large codebook ($>2000$ words) increased computational time and memory usage. Therefore, obtaining the most suitable values of BoW parameters for each specific data type produces better results while being more computationally efficient.

The generalization capability between the classification models for the discrimination of pain showed that the K-NN classifier performed slightly better than SVM on both channel sets. This can be explained due to the fact that the K-NN classifier was designed with a smaller codebook (500 words) compared to the SVM classifier (2000 to 2500 words). Thus, increasing the number of words amplifies the intraclass variability, which decreases classification accuracy and also increases computational time and memory usage. We also found that unlike SVM, K-NN was much more adept at classifying low-cold (threshold of cold) and high-heat (tolerance of heat) on both channel sets, which might be the reason why the accuracy rate for SVM was slightly lower compared to K-NN. The difference between these two modalities (low-cold and high-heat) of pain might also suggest that distinctive patterns from these two classes are more inconsistent among subjects and thus more difficult to separate in a high dimensional space using a linear kernel. As mentioned before, we also compared different values of NNs ($K=3$, 5, 7, 9), however, the performance dropped; the risk of including NNs that belong to different classes increases since the K-NN algorithm uses a majority voting scheme.

Multiple cortical areas were monitored using a 24-channel head probe, because of that, we were interested to know if the hemodynamic response in a specific RoI could provide better classification accuracy compared to the whole probe. For that reason, we selected the somatosensory area as RoI and trained the classifiers with data only from this cortical area using (eight channels). The reasoning behind this comparison was that the somatosensory area, in particular, the primary (S1) and secondary (S2) somatosensory cortex, has shown a fundamental role in pain processing.^8,12,18 The results showed that the accuracy of RoI-based classifiers was slightly lower than whole-probe classifiers (refer to Table 1). For instance, for the K-NN classifier, the difference between RoI-based and whole-probe approach was a loss of $\sim 0.55\%$ in accuracy, whereas for the SVM classifiers, the accuracy dropped 0.42%. A clear advantage of using an RoI-based approach is the need for less channels, which decreases the computational cost to build the BoW representation and the classification models (K-NN and SVM), since the RoI-based dataset is built with one third of the original dataset. However, the RoI-based approach did not perform better than the whole-probe approach in either case. This finding is consistent with Brown et al.,¹⁸ the authors found that using a whole-brain approach produced better accuracy to classify pain and no-pain signals using a linear SVM compared to an RoI-based approach. These findings confirm that despite the fact that the somatosensory region plays a fundamental role in pain processing, there are still more cortical regions than are implicated in sensory aspects of pain perception (e.g., insular cortex, anterior cingulate cortex, prefrontal cortex), which make nociceptive processing a large distributed brain network, referred as the “pain matrix.”⁴⁴ Therefore, it is also important to monitor different cortical regions that may be involved in the particular pain experience, the interaction between different cortical areas could tell more about the brain functioning as a system rather than looking at separate areas working independently.¹³

4.2.

Pain Assessment

The definition of pain is sometimes broad and controversial and is mainly regulated by the perceived intensity of the painful sensation. In this study, we differentiated the intensity of a painful stimulus as pain threshold (low pain) and pain tolerance (high pain). It is believed that all humans have a similar threshold for pain but different pain tolerance.⁴⁵ For instance, heat is perceived as noxious at thermal doses between 40°C and 45°C, same temperature at which tissue damage occurs, while the response to pain intensities differs from person to person and is influenced by factors, such as, gender, age, weight, etc. In addition, in some circumstances, the tolerance to pain of a person can change after a stressful experience (e.g., death of a loved one, self-injury). In a study by McCoy et al.,⁴⁶ it was found that participants with previous experience of self-injury had a significantly higher pain tolerance but the same pain threshold compared to normal controls; it was hypothesized that self-injury increases the desensitization to fear and pain.

An advantage of obtaining the intensity (low and high) of pain in a patient is to provide a personalized medical treatment. It is well reported that pain judgement is affected by factors, such as age, gender, weight, cultural background, etc. Medical treatments that work well for some patients may not be adequate for others due to physical differences, therefore tailoring pain control to the individual needs of each patient is not only important for nonverbal patients but also for the whole population. In many cases, medical practitioners prescribe analgesics based only on general information about medical treatments and not on the individual patient’s characteristics. For this reason, knowing the individual pain sensation will reduce the risks of harmful side effects due to overdose and cut medical treatment costs compared with traditional trial-and-error treatments.

Thermal testing was kept as uniform as possible among all the participants, however, possible differences in data acquisition could occur. In order to minimize any differences between observations and thermal acquisition, we utilized thermode-induced stimulation. In contrast with water-based methods, such as cold pressor test and the hot water immersion test,⁴⁷ the thermode induces heat or cold from the baseline of 32°C at a constant rate. One major difference between the two methods is the size of the area of application; while water-based methods require one hand to be submerged into cold/hot water, the thermode can be used on any area of the skin with a small surface ($9{\mathrm{cm}}^2$). Thermode-induced stimulus safely activates nociceptive peripheral afferents and induces a moderate to intense painful sensation. One important feature of thermode-based stimulation is its reproducibility, which reduces the intrasubject variability of the measurements and increases reliability of the thermal experiment. Thermode-induced stimulation is less time consuming, and it can be used routinely and repeatedly in patients, which makes it ideal for clinical use. In addition, thermode-based stimulation safely activates A-delta ($\mathrm{A}\delta$) and C fibre nociceptors, which are known to transmit sensory information to the spinal cord and the brain. Nociceptors are specialized sensory nerve endings that warn us of potential tissue damage by identifying harmful temperature, pressure, and chemicals; nociceptors convert these stimuli into electrical signals that are conveyed through the central nervous system to higher cortical regions in the brain. $\mathrm{A}\delta$ and C fibre nociceptors carry noxious sensory information and they have different functions while transmitting sensory information. Fine myelinated $\mathrm{A}\delta$ fibers respond to mechanical and cold stimuli, transporting rapid and sharp pain; unmyelinated C fibers respond to warm detection and heat pain, transmits slow, burning pain, responding to chemical, mechanical and thermal stimuli.²²

4.3.

Limitations

We acknowledge that one potential issue was the lack of control for any skin blood flow contributions and intracerebral hemodynamics to the fNIRS signals. Recent studies have highlighted the issue that fNIRS signals encompass not only hemodynamic fluctuations due to neurovascular coupling but also due to skin blood flow and task-related systemic activity of cortex.^23,28 On the one hand, extracerebral hemodynamics is originated at shallow tissue (e.g., scalp), the area in which NIR light travels before reaching the cerebral cortex, and thus, contributions to the observed fNIRS signals of HbO and Hb are expected. Therefore, it is important to account for probable influence of scalp blood flow, which may or may not be temporally consistent with cortical signals. In our study, despite possible contributions of extracerebral confounders, the observed activation was measured in areas known to be involved in pain processing, which suggests that the fNIRS signals reflect mainly localized cortical vascular dynamics. On the other hand, intracerebral hemodynamics is caused by systemic physiological interference that may resemble true task-related cortical activities. These physiological variables are mainly associated with cardiac pulsations, respiration, arterial blood pressure, and activity in the autonomic nervous system.²³ In our case, the peak amplitudes for the HbO response to thermal stimuli were substantially distinguishable relative to the baseline signal at resting state. In both cases, we tried to reduce these confounders (extra- and intracerebral) by (1) using PCA on the baseline data to identify those PC that account for most of the variance while in resting state, and then reduce these PCs (potentially both, extra- and intracerebral confounders) in the stimulus data. (2) Using wavelet analysis in frequency bands that are only related to the expected period (Hz) of response and therefore avoiding those bands related to heartbeat ($\sim 1\mathrm{Hz}$), respiration ($\sim 0.3\mathrm{Hz}$), or blood pressure ($\sim 0.1\mathrm{Hz}$); however, very-low frequency oscillations (VLFO) associated with various brain stem nuclei in the sympathetic tone of cerebral arterioles (0.02 to 0.08 Hz)⁴⁸ might overlap with the content of the stimulus response, since the period of our block-design stimulation ($\sim 30\mathrm{s}$) corresponds to this band. Finally, using both methods might not be enough to reduce such sources of hemodynamic interference and other methods, such as using short-separation channels to reveal scalp blood flow, or employing additional instruments to measure physiological signals would be desirable to regress out of these global signals.

4.4.

Future Work and Final Remarks

The results of this study expand earlier findings from previous research of pain assessment using neuroimaging methods.^18,19,42 Direct comparisons are difficult because of the use of different experimental conditions, neuroimaging technique (e.g., fMRI), sampled population, data representation (using BoW), validation method, and classification models. Therefore, the contributions of this study can be summarized as follows: (1) presenting the classification of two types of thermal stimuli (cold and heat) and the discrimination between corresponding low (threshold) and high (tolerance) pain for each thermal stimulus, (2) showing that the pain tolerance and pain threshold test from the QST could be used to obtain pain information from nonverbal patients, (3) comparing two classification models (K-NN and SVM) to classify four different modalities of pain, and (4) investigating whether feature extraction from the somatosensory area could produce better classification accuracy compared to whole-probe approach.

The development of a bedside monitor for the diagnosis of pain is yet to be fully established. To expand such method into clinical settings requires validation on new datasets (generalization), include more participants of different populations (e.g., age, background), use multiple types of noxious stimulations (e.g., mechanical, chemical), stimulate different parts of the body (robust to testing location), and include models to study chronic pain disorders (e.g., back pain, neuropathic pain, fibromyalgia). In addition, the use of more learning models should be evaluated, such as decision trees, linear discriminant analysis (LDA), neural networks (NN), or higher-order SVM kernels (e.g., radial basis function, polynomial).

fNIRS has demonstrated to be a method that has potential for the assessment of pain. fNIRS is a technique capable of identifying cortical hemodynamic changes in response to chemical, temperature, and pressure noxious stimuli. In addition, our study demonstrates that the use of fNIRS in combination with machine learning techniques is a powerful tool for the assessment of pain in experimental settings. fNIRS possesses advantages over PET or fMRI to be used in more realistic clinical settings, e.g., less exposure to ionizing radiation, safe to use over long periods and many times, less expensive, easy to use, and small size. Finally, our findings advance knowledge in pain assessment using neuroimaging as a method of diagnosis and represent a step closer to developing a physiologically based diagnosis of human pain that would benefit vulnerable population who cannot self-report pain.