Pain Detection with fNIRS-Measured Brain Signals: A Personalized Machine Learning Approach Using the Wavelet Transform and Bayesian Hierarchical Modeling with Dirichlet Process Priors

TL;DR

This study develops a personalized machine learning method using wavelet features and Bayesian hierarchical modeling with Dirichlet process priors to detect pain from fNIRS signals, aiming for objective pain assessment.

Contribution

It introduces a novel personalized approach combining wavelet transform and Bayesian hierarchical modeling for pain detection using fNIRS data.

Findings

Effective pain detection from prefrontal cortex fNIRS signals.

Personalized modeling improves detection accuracy.

Supports use of fNIRS as an objective pain assessment tool.

Abstract

Currently self-report pain ratings are the gold standard in clinical pain assessment. However, the development of objective automatic measures of pain could substantially aid pain diagnosis and therapy. Recent neuroimaging studies have shown the potential of functional near-infrared spectroscopy (fNIRS) for pain detection. This is a brain-imaging technique that provides non-invasive, long-term measurements of cortical hemoglobin concentration changes. In this study, we focused on fNIRS signals acquired exclusively from the prefrontal cortex, which can be accessed unobtrusively, and derived an algorithm for the detection of the presence of pain using Bayesian hierarchical modelling with wavelet features. This approach allows personalization of the inference process by accounting for inter-participant variability in pain responses. Our work highlights the importance of adopting a personalized approach and supports the use of fNIRS for pain assessment.

Tables

Table 1. TABLE I: Classification performance of single-task models

Model	Accuracy	Pr.	R.	F1
Logistic regression (L1)	0.58	0.59	0.55	0.55
Logistic regression (L2)	0.60	0.61	0.57	0.58
SVM (linear kernel)	0.57	0.58	0.55	0.55
SVM (rbf kernel)	0.69	0.72	0.63	0.67

Table 2. TABLE II: Classification performance of the HBLR model

K	Accuracy	Pr.	R.	F1
2	0.66	0.53	0.53	0.53
3	0.73	0.67	0.53	0.59
4	0.81	0.75	0.71	0.73
5	0.75	0.67	0.59	0.63

Equations

$${\Psi }_{w_o}(t)=c_{w_o}{\pi }^{-\frac{1}{4}}{e}^{-\frac{1}{2}{t}^2}\left(e^{i{w}_{o}t}-e^{-\frac{1}{2}{w}_o^2}\right)$$

$$T(a,b)=\frac{1}{\sqrt{a}}{\int }_{-\infty }^{+\infty }x(t){\psi }^{\ast }\left(\frac{t-b}{a}\right)dt$$

$$p\left(y_{m,i}|{\mathbf{w}}_m,{\mathbf{x}}_{m,i}\right)=\sigma {\left({\mathbf{w}}_m^T{\mathbf{x}}_{m,i}\right)}^{y_{m,i}}{\left[1-\sigma \left({\mathbf{w}}_m^T{\mathbf{x}}_{m,i}\right)\right]}^{1-y_{m,i}}$$

$$w_m|G$$

$$G$$

$$p\left(\alpha |{\tau }_{10},{\tau }_{20}\right)=\frac{{\tau }_{20}^{{\tau }_{10}}}{\Gamma \left({\tau }_{10}\right)}{\alpha }^{{\tau }_{10}-1}\exp \left(-{\tau }_{20}\alpha \right)$$

$$p\left(\mathbf{Z}|{\left\{{\mathcal{D}}_m\right\}}_{m=1}^M,\mathbf{\Phi }\right)=\frac{p\left({\left\{{\mathcal{D}}_m\right\}}_{m=1}^M|\mathbf{Z},\mathbf{\Phi }\right)p(\mathbf{Z}|\mathbf{\Phi })}{p\left({\left\{{\mathcal{D}}_m\right\}}_{m=1}^M|\Phi \right)}$$

$$q(\mathbf{Z})=\left[{\prod }_{m=1}^M{q}_{c_m}\left(c_m\right)\right]\cdot \left[{\prod }_{k=1}^K{q}_{v_k}\left(v_k\right)\right]\cdot q_{\alpha }(\alpha )\cdot [\left({\prod }_{k=1}^K{q}_{w_k^{\ast }}\left(w_k^{\ast }\right)\right]$$

$$p(y_{m,\ast }=1|x_{m,\ast },\Phi ,$$

$$=$$

$$\approx$$

Full text

I Introduction

Functional near infrared spectroscopy (fNIRS) is an emerging brain-imaging technique that provides non-invasive, long-term measurements of cortical hemoglobin concentration changes [1]. With fNIRS, near infrared light is projected onto the surface of the head, and is mainly absorbed by two types of hemoglobin during its propagation within the cerebral cortex, i.e. oxygenated hemoglobin (HbO) and deoxygenated hemoglobin (HbR). By continuously monitoring attenuation of light, fNIRS is able to reconstruct the cortical concentration changes of HbO and HbR, from which the local neuronal activity can be inferred. Compared with other imaging modalities such as functional magnetic resonance imaging (fMRI), fNIRS shares a similar physiological basis but has significant advantages in cost, robustness and portability, making it suitable to be used especially at bedside or in complex clinical settings.

One novel application of fNIRS is to provide objective and robust assessment of pain. Pain is a subjective and complex experience that is processed within the human nervous system [2]. While the current methods of assessing pain mostly rely on the participants’ self-report (the “gold standard”) [3] or physiological signals (e.g. heart rate, blood pressure or skin conductance) [4, 5, 6, 7, 8, 9], establishing an objective, robustness marker of pain in the central nervous system with a portable brain-imaging device would be beneficial for participants who cannot effectively communicate. In particular, patients undergoing surgery are unconscious during actual tissue damage. Without proper detection and control, repetitive pain signaling during surgery is postulated to lead to severe pain in the postoperative period and also potential initiation of chronic neuropathic pain [10, 11].

Several previous studies employing fNIRS have revealed specific patterns in the cortical response to evoked pain stimuli from healthy volunteers, which were mainly characterized by activations (i.e. increases in HbO concentration and decreases in HbR concentrations, inferring elevated local neuronal activity) in the sensorimotor cortex (SMC) as well as deactivations (i.e. decreases in HbO concentration and increases in HbR concentration, inferring inhibited local neuronal activity) in the medial prefrontal cortex (mPFC) [12, 13]. In our recent work exploring the possibility of combining fNIRS and machine learning to detect pain, we extracted features from both the SMC and the mPFC signals, and applied multi-task multiple kernel machine learning to account for the inter-participant variability in the brain response to pain [14]. While this approach yielded moderate detection accuracy ($\sim$80%), it is difficult to implement in the operating room, mainly because of the challenges of putting fNIRS optodes over the SMC in surgical patients who are usually kept in a supine position. Therefore, in this study, we decided to focus only on the prefrontal signals.

Previous work on pain detection has highlighted the importance of adopting a personalized machine learning approach to account for person-specific differences in pain responses [7, 8, 15, 14]. Here, we used multi-task learning (MTL) to learn personalized classifiers for the presence of pain while leveraging data from the entire population [16]. Specifically. we employed a non-parametric Bayesian hierarchical model to learn personalized logistic regression classifiers [17] from time-frequency features derived using the continuous wavelet transform [18]. Our fNIRS-based machine learning approach was evaluated on a dataset containing pain responses to electrical noxious stimuli as well as data from baseline states.

The main contributions of this work are: (1) we present an approach to identify pain responses using fNIRS, (2) only fNIRS channels from the prefrontal cortex were used, (3) we show that by personalizing the machine learning models using multi-task learning the system is able to provide differentiation of pain state with better performance. This approach could have great promise for pain assessment in noncommunicative patients as well as multiple clinical populations, and form the basis of a more objective measure of pain than self-report.

II Data

II-A Data acquisition

This study was approved by the institutional review boards of Massachusetts General Hospital and Boston Children’s Hospital. The datasets of a total of 43 healthy participants (age range: 26.8 $\pm$ 5.6, 20 to 39 years old) were included in this study. Previous work has shown sex differences in pain-related cerebral processing and blood oxygenation signals [19, 20]. Hence, in this work we decided to focus on a single gender to control the number of varying factors in the data. Because pain-related brain activation in females is also be influenced by menstrual cycles [21], only male participants were enrolled. Written consent was obtained from each of the participating participants prior to the study procedures.

A TechEn CW7 continuous wave fNIRS system (TechEn Inc., Massachusetts, USA) was used for data acquisition. Three light emitters and 4 light detectors were mounted to cover the prefrontal areas of the brain, forming 8 emitter-detector pairs (called fNIRS channels; see Fig.1) with inter-optode distance of 3 cm to enable detection sensitivity of cortical hemodynamic activities. Two wavelengths of near infrared light were used: 690 nm and 830 nm. One intrinsic limitation of fNIRS is that the recorded data contains not only the hemodynamic signals from the cerebral cortex but also the contributions from several extracerebral layers (e.g. skin, scalp and skull) because near infrared light has to penetrate those layers before and after its travel within the cortex. Those extracerebral signals are mostly related to the non-neuronal physiological effects of the body including heartbeat, respiration and blood pressure. To correct for such contamination, we installed a short distance detector 8 mm away from each of the light emitters to just sample the hemodynamic changes from those extracerebral layers. The recorded short distance data were later used in the general linear model (GLM) analysis as regressors (see below).

During the data acquisition period, each participant underwent the following three scanning sessions supervised by study personnel: a resting-state session, a tactile brush session and an electrical stimulation session. During the resting state, the participant was asked to simply sit comfortably in a chair and rest for 6 min. In the tactile brush session, a study personnel delivered 12 non-painful brush stimuli by gently brushing the dorsum of the participant’s left hand. Each brush stimulus lasted for 5 sec and was followed by a resting period of 25 sec. In the electrical stimulation session, two levels of 5 Hz electrical pulses were applied by the study personnel to the participant’s left thumb with a neurometer (Neurotron, Maryland, USA): an innocuous level (score: 3/10) for which the participant described that he was strongly aware but did not feel any pain, as a noxious level (score: 7/10) for which the participant reported much pain but was able to tolerate without breath holding, sweating or any retreat actions. The intensities of the two electrical stimulation levels were determined with a few pre-tests before the actual scanning session. Similar to the tactile brush stimulus, each electrical stimulus lasted 5 sec and was separated from one another with a 25 sec resting period. The entire electrical stimulation sequence contained a total of 6 innocuous stimuli and 6 noxious stimuli. The order of the innocuous and noxious stimuli was randomized for each participant to accommodate for the effect of pain anticipation.

II-B Data preprocessing

The acquired fNIRS data were first processed with an open-source toolbox HomER2 [22] based on Matlab (Mathworks, Massachusetts, USA). The measured light intensity of the channels (see Fig.1) was then converted to optical density changes by taking the logarithm of the signal. Motion artifacts were then detected and corrected with the hybrid Spline Interpolation and Savitzky-Golay filtering method [23]. A temporal low-pass filter with a cutoff frequency of 0.5 Hz was applied to remove high frequency oscillations in the data that presumably did not have a neuronal basis (e.g. heart beat whose frequency is around 1 Hz). The processed optical density timecourses were then transformed to the concentration changes of HbO and HbR using the modified Beer-Lambert Law [24, 25, 26] with a partial pathlength factor of 6. The concentration changes of total hemoglobin (HbT) were obtained by a direct summation of the HbO and HbR changes (HbT = HbO + HbR). For each of the tactile brush and electrical stimulation sessions, we then reconstructed the hemodynamic response function (HRF) to tactile brush or to noxious and innocuous electrical stimuli from 2 sec before the onset of the stimulus to 20s after the onset for each of the three types of hemoglobin using a GLM [27]. In the model, the HRFs to stimuli were approximated with consecutive Gaussian basis functions with width and temporal inter-spacing to be both 1 sec. For each normal fNIRS channel, the short separation measurement that had the highest correlation was also included in the GLM as a regressor to filter out the physiological signal from extracerebral layers. Polynomial terms up to the 3rd order were also added to the model to remove the low frequency drifts from the signal.

II-C Feature extraction

Following [14], we extracted windows of duration 20 seconds from the HbO signals. These windows were either ”pain” or ”no-pain” windows. Pain windows were extracted starting immediately after the onset of the electrical noxious stimulus (score: 7/10) and were assigned label $y=1$. No-pain windows were randomly sampled from the pre-stimulation baseline (resting) recording and were assigned label $y=0$.

Feature extraction was performed using the discretized continuous wavelet transform [18], a powerful time-frequency analysis technique that provides a description of the power spectrum of the signal in terms of both time and frequency domains simultaneously. Hence, it remedies the drawbacks of traditional time-domain feature extraction (inability to characterize the frequency components of the fNIRS signals) and frequency-domain feature extraction (temporal dynamics are lost).

Wavelets are functions with zero mean that are both time and frequency localized. In this work, we use the Morlet Wavelet, which is given by:

[TABLE]

where $c_{w_{o}}=\left(1+e^{-w_{o}^{2}}-2e^{-\frac{3}{4}w_{o}^{2}}\right)^{-\frac{1}{2}}$ is a normalization constant, $w_{0}$ is the wavelet central frequency and $t$ is time [28]. The wavelet transform is defined as follows:

[TABLE]

where $x(t)$ is the signal analyzed, $a$ is the dilation of the wavelet, $b$ is the location of the wavelet and $\psi^{*}(t)$ is the complex conjugate of the wavelet function $\psi(t)$. This continuous wavelet transform was discretized using ten voices per octave.

Following the wavelet transformation, we extracted the scalogram (see Fig.2 for an example), which represents the strength or energy of each coefficient, using the $L_{2}$ norm. From the scalogram, we extracted features from two frequency bans: very-low frequency oscillations (VLFO, 0.01-0.08 Hz) and low frequency oscillations (LFO, 0.08-0.15 Hz). Specifically, we extracted the following features in the wavelet domain: (1,2,3) the mean, maximum and standard deviation of $|T(a,b)|$ in $a=[0.01,0.08]$, (4,5,6) the mean, maximum and standard deviation of $|T(a,b)|$ in $a=[0.08,0.15]$, (7,8) the location of the maximum and slope of $\int_{a=0.01}^{0.08}|T(a,b)|da/\int_{a=0.01}^{0.08}da$, (9,10) the location of the maximum and slope of $\int_{a=0.08}^{0.15}|T(a,b)|da/\int_{a=0.08}^{0.15}da$. This was done for each of the 8 fNIRS channels depicted in Fig.1, resulting in a feature vector of dimension $D=80$. All features were then normalized prior to being used in our ML models described in Sec.III.

III Personalized machine learning model

In this work, we consider the binary classification task of detecting the presence of pain using multi-task learning. Traditional single-task learning (STL) refers to the approach of learning a classifier for a single classification task using the corresponding dataset $\mathcal{D}=\{(\bm{x}_{i},y_{i})\}_{i=1}^{N}$, where $\bm{x}_{i}\in\mathbb{R}^{D}$ is a $D$-dimensional feature vector, $y_{i}\in\{0,1\}$ is a binary label and $N$ is the number of samples in the dataset. Many classification tasks can be viewed as consisting of multiple correlated subtasks, such as pain detection in a particular patient subpopulation that responds to noxious stimuli in a characteristic way, similar to yet different from other patient subpopulations [14, 7, 8, 15]. Therefore, pooling all subtasks and treating them as a single task may not be appropriate. Furthermore, the traditional STL approach of isolating each subtask and learning the corresponding classifier independently may not be appropriate neither, as it does not exploit the potential information that may be acquired from other correlated classification subtasks.

Multi-task learning (MTL) is a type of inductive transfer learning in which multiple classifiers are learned simultaneously for several related tasks using a shared representation [16]. This is distinct from traditional single-task learning in several aspects: (1) multiple tasks are learned in parallel, which may represent e.g. individuals or groups of individuals that share certain characteristics, (2) these tasks are not identical, so that pooling them together into a single task is not appropriate, (3) because tasks are correlated, what is learned form one task may be transferable to another, (4) by learning tasks in parallel, tasks can exploit the potential information one may acquire from other classification tasks, hence leveraging the limited amount of training data available for each task, to the benefit of all. Furthermore, exploiting data from related tasks leads to better generalization of the resulting models.

In this work, following previous work on multi-task learning for predicting mood, stress, and health [29], we used hierarchical bayesian logistic regression (HBLR). This a nonparametric hierarchical Bayesian MTL model with a common prior drawn from the Dirichlet process (DP) that learns logistic regression classifiers. This model, introduced in [17], results in personalized logistic regression classifiers that account for individual characteristics in brain hemodynamic responses to pain, while learning from the entire population. Moreover, our approach induces a soft clustering of the tasks without any a priori knowledge or meta-information. This may be used to uncover novel pain phenotypes.

Let’s consider $M$ tasks with corresponding datasets $\mathcal{D}_{m}=\{(\bm{x}_{m,i},y_{m,i})\}_{i=1}^{N_{m}}$, where $m\in\{1,...,M\}$ is the task and $N_{m}$ is the number of instances for task $m$. Following [17, 29], we model the conditional distribution of $y_{m,i}$ given $\bm{x}_{m,i}$ via logistic regression:

[TABLE]

where $\sigma(x)=\frac{1}{1+\exp{(-x)}}$ and $\bm{w}_{m}$ represents the classifier weights for task $m$. The goal is to learn $\{\bm{w}_{m}\}_{m=1}^{M}$ jointly, sharing information between tasks. To do so, the hierarchy consists of a bottom layer with task-specific parameters, and a top layer in which tasks are connected together via a common prior. Given the prior, individual models are learned independently. However, the common prior is also learned during training, resulting in information transfer between tasks.

As in [17], the common prior $G$ on the task-specific model parameters in the proposed nonparametric hierarchical Bayesian model is drawn from the Dirichlet process (DP):

[TABLE]

where $\alpha\sim Ga\left(\tau_{10},\tau_{20}\right)$ is the positive innovation parameter drawn from a Gamma distribution and $G_{0}\sim N_{d}(\bm{\mu}=\bm{0},\bm{\Sigma}=\sigma\bm{I})$ is $d$-dimensional multivariate normal base distribution. The DP prior is employed due to its implicit non-parametric soft-clustering mechanism, which enables the model to automatically identify the similarities between the various tasks and adjust the complexity of the model, that is, the number of task clusters. This also provides valuable insights into tasks that exhibit similar behaviours. Because the clustering is soft, tasks need not be assigned discretely to single clusters. Instead, tasks can belong to many clusters in varying degrees. The innovation parameter $\alpha>0$ controls the probability of creating a new cluster, with larger $\alpha$ yielding more clusters. When $\alpha\rightarrow\infty$ there is a cluster for each task, whereas small values of $\alpha$ lead to only a few distinct clusters. In this work, $\alpha$ has the following probability distribution function in the shape-rate parametrization:

[TABLE]

where $\Gamma(\cdot)$ is the Gamma function. Integrating $\alpha$ over a diffuse hyper-prior, as opposed to defining it directly, increases the roboustness of the algorithm [17].

In Bayesian modeling, we are interested in the posterior distribution of the latent variables given the training data and hyperparamters:

[TABLE]

where $\bm{Z}=\{\{\bm{c}_{m}\}_{m=1}^{M},\{v_{k}\}_{k=1}^{\infty},\alpha,\{\bm{w}_{k}^{*}\}_{k=1}^{\infty}\}$ denotes the collection of latent variables, $\bm{\Phi}=\{\tau_{10},\tau_{20},\bm{\mu},\bm{\Sigma}\}$ denotes the collection of given parameters and hyper-parameters and $\bm{c}_{m}$ is an all-zero vector except that the $k$-th entry is equal to one if task $m$ belongs to cluster $k$. To approximate $p\left(\mathbf{Z}|\left\{\mathcal{D}_{m}\right\}_{m=1}^{M},\Phi\right)$, whose computation is intractable, we used mean-field variational bayesian (VB) inference [30]. This approximates the true posterior by a variational distribution $q(\bm{Z})$ and converts computation of posteriors into an optimization problem. Following [31], we adopted a truncated stick-breaking representation for the variational distribution. By setting the truncation level equal to the number of desired clusters $K$, we can control the number of resulting clusters and the complexity of the HBLR model. The factorized variational distribution is then specified as:

[TABLE]

where $c_{m}\sim M_{K}\left(1;\phi_{m,1},\ldots,\phi_{m,K}\right),m=1,\ldots,M$, $v_{k}\sim Be\left(\varphi_{1,k},\varphi_{2,k}\right),k=1,\ldots,K-1$ and $w_{k}^{*}\sim N_{d}\left(\theta_{k},\Gamma_{k}\right),k=1,\ldots,K$ (see [17] for details).

To learn all the parameters, we use a coordinate ascent algorithm [17] that uses the mean-field approach [30]. All hype-parameters were initialized following [29], and re-estimated iteratively until convergence.

The prediction function for a new test sample $x_{m,*}$ becomes:

[TABLE]

where we use the approximate form of the integral [32] and $\sigma(x)=1/(1+e^{-x})$ is the sigmoid function.

IV Results

After processing the fNIRS signals as described in Sec.II-B, windows of duration 20 seconds were extracted from the HbO signals. The choice of signal modality (that is, HbO) and window size was informed by previous work done on pain detection from fNIRS [33, 14]. From these windows, we extracted the $D=80$ features described in Sec.II-C from the prefrontal fNIRS channels depicted in Fig.1. All features were normalized, and the dataset was balanced by downsampling the over-represented class (no pain; $y=0$), resulting in $N=510$ instances from $M=43$ tasks with both classes equally represented.

First, we evaluated traditional single-task machine learning models using 10-fold cross-validation. Namely, we used logistic regression with $L_{1}$ and $L_{2}$ regularization [34], and support vector machines with linear and radial basis function (RBF) kernels. The later resulted in the highest accuracy (0.69) among these models, as shown in Table I.

Following the single-task models, we shifted to the multi-task setting and evaluated the multi-task HBLR model with different values of $K$, which corresponds to different truncation levels or number of clusters. To do so, for each of the 43 tasks we randomly split the data into training (90%) and test partitions (10%) using 10-fold cross-validation. The results are shown in Table II and indicate improved performance of the HBLR model with respect to single-task models. The best performance was achieved with $K=4$, and significantly outperforms the single-task models. We also evaluated $K>5$ (not shown in Table II), but the additional number of clusters did not result in improved performance. For every value of $K$, model hyperparameters were optimized using 10-fold cross-validation. The best performing model used $\tau_{10}=0.01$, $\tau_{20}=0.1$.

The HBLR model performs a soft clustering of the $M=43$ tasks into $K$ clusters, such that a given task may belong to multiple clusters with varying degrees of membership. This is visualized in Fig.4 for the best performing model with $K=4$. While the analysis of the clustering results may be used to uncover novel pain phenotypes, this was hampered in this work by the limited dataset size, low heterogeneity in the participant pool, as well as lack of meta-information.

V Discussion and Conclusion

This study investigated the use of functional near-infrared spectroscopy (fNIRS) for the detection of evoked pain, using fNIRS signals derived from optodes placed exclusively on the frontal brain cortex. Multiple studies have pointed towards an integrated role of the prefrontal cortex in the high-level process of pain perception [14, 35]. Whereas previous studies (e.g. [14]) have included other brain regions such as the sensorimotor cortex, focusing only on prefrontal cortical area has many advantages including its ease of access (especially in surgical settings where the patient is kept in a supine position) and the ability to avoid hair contamination in the recorded fNIRS signals. Moreover, the need of only sampling from the prefrontal cortex greatly reduces the required number of fNIRS optodes and optical fibres to be mounted on the participant’s head, leading to an easier installation process and more patient comfort. Hence, this work investigated the use of frontal fNIRS channels exclusively. One central challenge in establishing objective and reliable assessment of pain with brain signals is the notable inter-subject variability in pain perception and responses, which hampers the ability of machine learning models to generalize across people. Previous work has shown improved pain assessment performances with personalized machine learning models [7, 8, 15, 14]. Hence, this work employed multi-task machine learning to provide personalized analysis.

We evaluated our model on a dataset containing fNIRS responses to experimental electrical pain. Our results confirmed previous work [14] that showed that fNIRS-measured brain signals and machine learning may be used to detect the presence of pain. Using exclusively frontal channels, we showed that high pain state discrimination accuracy may be achieved. Furthermore, our results confirmed the importance of adopting a multi-task approach to personalize the assessment.

This work has several limitations. The strength of the model is limited by the size of our dataset (43 participants) and lack of heterogeneity in the participant pool of healthy males. Moreover, only the detection of experimental evoked electrical pain at a single pain intensity was evaluated. Future work should extend our model to females, other participant phenotypes, different pain stimuli, varying pain intensities, and other types of pain (e.g. ongoing or chronic pain). Also, the model was evaluated using pain windows aligned with the onset of the stimulus. Future work should evaluate the model continuously (see [36] for an example). Despite the limitations, this novel fNIRS personalized machine-learning approach to pain detection showed great promise in its accuracy, using only prefrontal electrodes. Thus, this work advances progress toward objective pain recognition, providing accurate and practical measures of pain that work even if patients lack alertness or the ability to articulate their experience. It may also help advance the development of automatic analgesia administration systems for hospital settings that currently rely on subjective self-reported pain measures [37].