arXiv:1711.04036v1 [] 10 Nov 2017

Physiological and behavioral profiling for nociceptive
pain estimation using personalized multitask learning

Daniel Lopez-Martinez1,2 , Ognjen (Oggi) Rudovic2 , Rosalind Picard2
1
Harvard-MIT Health Sciences and Technology, dlmocdm@mit.edu
2
Affective Computing group, MIT Media Lab, Massachusetts Institute of Technology

Abstract
Pain is a subjective experience commonly measured through patient’s self report.
While there exist numerous situations in which automatic pain estimation methods
may be preferred, inter-subject variability in physiological and behavioral pain
responses has hindered the development of such methods. In this work, we address
this problem by introducing a novel personalized multitask machine learning
method for pain estimation based on individual physiological and behavioral pain
response profiles, and show its advantages in a dataset containing multimodal
responses to nociceptive heat pain.

1

Introduction

Pain is an unpleasant sensory and emotional experience associated with actual or potential tissue
damage with sensory, emotional, cognitive and social components [1]. In the clinical and research
settings, pain intensity is measured using patient’s self-reported pain rating scales such as the visual
analog scale (VAS) [2]. Unfortunately, these self-report measures only work when the subject is
sufficiently alert and cooperative, and hence they lack utility in multiple situations (e.g. during
drowsiness) and patient populations (e.g. patients with dementia or paralysis).
To circumvent these limitations, automatic methods for pain estimation based on physiological
autonomic signals [3, 4] and/or facial expressions [5, 6] have been proposed. However, inter-subject
variability in pain responses has limited the ability for the automated methods to generalize across
people. For example, autonomic responses captured in signals such as heart rate and skin conductance
have been found to be correlated only moderately with self-reported pain [7]. A large part of the
variance may be explained by inter-individual variability in autonomic reactivity, independent from
stimulus intensity, and also by inter-individual variability in brain activity within structures involved
in regulation of nociceptive autonomic responses [7]. Similarly, facial responses [8] also vary across
individuals [9, 10]: Some people show strong facial expressions for very low pain intensities, while
others show little or no expressiveness. Therefore, it is important to account for individual differences.
While several recent works have shown the advantages of personalization for pain estimation, both
from physiological signals [3] and from face images [6, 5], none of these approaches explored the
effect on pain estimation performance of clustering subjects into different clusters or profiles. Hence,
in this work we investigate the grouping of subjects into different profiles based on their unique
multimodal physiological and behavioral (facial expression) responses to pain. These profiles are
then used to define the structure of a multi-task neural network (MT-NN), where each task in the
MT-NN corresponds to a distinct profile. The advantages of this work are two-fold: (i) we show that
the proposed multimodal approach achieves better performance than single-modality approaches, and
(ii) we also show that by accounting for the different profiles using multi-task machine learning we
achieve further improvements in pain estimation compared to single-task (population level) models.
31st Conference on Neural Information Processing Systems (NIPS 2017), Long Beach, CA, USA.

2

Dataset, Methods and Experiments

We used the publicly available BioVid Heat Pain database [11], which contains 87 participants
with balanced gender ratio. The experimental setup consisted of a thermode that was used for pain
elicitation on the right arm. Before the data recording started, each subject’s individual pain threshold
(the temperature for which the participant’s sensing changes to pain) and tolerance threshold (the
temperature at which the pain becomes unacceptable) were determined. These thresholds were used
as the temperatures for the lowest and highest pain levels together with two additional intermediate
levels, thus obtaining four equally distributed pain levels (P = {1, 2, 3, 4}). For each subject, pain
stimulation was applied 20 times for each of the 4 calibrated temperatures, for 4 seconds followed
by a recovery phase (8-12 seconds, randomized). Each recording lasted 25 min on average. The
following signals were recorded: (1) skin conductance, (2) electrocardiogram, and (3) face videos.
2.1

Signal processing and feature extraction

Our pain estimation approach is posed as a continuous regression problem. Hence, for each recording
we sample windows of duration 6 seconds with a step size of 0.5 seconds, and assign the label
corresponding to the heat applied at the beginning of the window. To build the profiles, we used
windows of duration 8 seconds, as described in Sec. 2.2. In what follows, we describe the features
extracted for each window.
Skin conductance (SC): We used nonnegative deconvolution [12] to decompose the SC signal into
its tonic and phasic components, and to extract the phasic driver, a correlate of the activity of the
sudomotor nerves [13], such that SC = SCtonic + SCphasic = SCtonic + Driverphasic ∗ IRF, where IRF
is the impulse response function of the SC response (SCR). Then, for each window, we extracted
the following features: (1) the number of SCRs with onset in the window, (2) the sum of amplitudes
of all reconvolved SCRs (SCphasic ) with onset in the window, (3) average phasic driver activity, (4)
maximum phasic driver activity, (5) integrated phasic driver activity, and (6) mean tonic activity.
Electrocardiogram (ECG): From the ECG, we extracted the R peaks [14] and subsequently
calculated the inter-beat intervals (IBIs). Based on previous work [15], the following features were
extracted for each window: (1) the mean of the IBIs, (2) the root mean square of the successive
differences, (3) the mean of the standard deviations of the IBIs, and (4) the slope of the linear
regression of IBIs in its time series.
Face: We used OpenFace [16], an open source facial behaviour analysis toolkit, to extract for each
frame the following features from the video sequences: locations of 68 facial landmarks, 3D head
pose and eye gaze direction [17], and 17 facial action units (AUs) [18]. From the facial landmark
locations, we extracted a set of geometric-based features. To do so, we first registered all landmarks
using the affine transform with 4 reference points: the two eye centers, the nose center and the
mouth center. Then, as in [19], at the frame level, we calculated the euclidean distance between
each of the 41 facial landmarks (we excluded the face contour and eyebrows) and the center of
gravity of the facial landmark set. Therefore, each frame was represented by a 41D vector. Then, at
the window level, we calculated 4 statistical features for each distance (mean, standard deviation,
max, min), hence obtaining 164D geometric-based features for each window. We also calculated the
same statistical descriptors for the eye gaze coordinates (24D) head pose (24D), and intensities of 17
actions units (AUs) detected by OpenFace, which are given on a scale from 0-5.
2.2

Subject profiling with spectral clustering

To build the profiles, we extracted windows corresponding to the first 48 heat stimuli. Specifically,
we used windows of duration of 8 sec., starting immediately after the onset of the stimulus. To profile
the subjects based on their physiology, we calculated the SC and ECG features described in Sec. 2.1.
For facial behavioral profiles, we computed the facial expressiveness of the subjects, defined as the
amount of variability in the landmark coordinates within the 8 sec. windows. From each sequence,
we calculated average
coordinates of each

 of the registered facial landmarks li = (xi , yi ), such that
¯li = (x̄i , ȳi ) = PF xf,i , PF yf,i , where F is the number of frames in each 25 min sequence,
f =1
f =1
and i ∈ [1, 68], provided by OpenFace [16]. Then, for each landmark i and each frame f = 1, . . . , N ,
in a given window, we calculated the distance from the landmark coordinates to the mean coordinates,
2

Figure 1: Subject profiling using spectral clustering. The i, j elements in the (clustered) similarity
matrices represent the distance from person i to person j in their physiological and behavioral features.
The more yellow the matrix element, the closer the corresponding subjects are.
# clusters
NC
c=2
c=3

Cluster
1
1
2
1
2
3

Size
85
40
45
23
35
27

Males
50.59%
55.00%
46.67%
56.52%
48.57%
48.15%

Females
49.41%
45.00%
53.33%
43.48%
51.43%
51.85%

Age (mean(std))
41.20(14.36)
44.53(13.11)∗
38.24(14.78)∗
47.61(11.24)∗
38.03(14.80)
39.85(14.42)

Table 1: Cluster statistics for the clusters shown in Fig. 1. Figure 2: Example of face tracked with
For each cluster, statistical significance (p-value ≤ 0.5) in OpenFace [16] showing landmarks,
terms of age composition with respect to the other clusters head pose and eye gaze.
is indicated with ∗.
i.e., di,f = ||li,f − ¯li ||. Finally, for each window, the facial expressiveness
scoreis obtained as the
PN P68
1
level of variability in facial landmarks, computed as lexp = N f =1
i=1 di,f .
To obtain the subject descriptors, for each subject we computed the mean value of each of the above
11 features, within windows of each of the 4 pain levels separately. Hence, each subject s was
represented by a 44D vector ps = [ps,1 , ..., ps,44 ]. We normalized these vectors to sum to 1 across
the pain levels per subject (to compensate for the scale differences among subjects), and used them to
cluster subjects into different profiles. To this end, we used normalized spectral clustering [20, 21].
Namely, let p̂s be the normalized descriptor of subject s. First, we construct a fully connected
similarity graph and the corresponding weighted adjacency matrix W . We used the fully-conected
graph, with edge weights wij =Ki,j = K(pi , pj ), and the radial basis function (RBF) kernel:
K(pi , pj ) = exp −γ||pi − pj ||2 with γ = 0.18, as the similarity measure. Then, we build the
degree matrix D as the diagonal matrix with degrees d1 , ..., dn on the diagonal, where di is given by
PM
di = j=1 wij where M is the number of subjects in our dataset. Next, we compute the normalized
graph Laplacian L = I − D−1 W and calculate the first c eigenvectors u1 , ..., uc of L, where c is the
desired number of clusters. Let U ∈ RM ×c be the matrix containing the vectors u1 , ..., uc as columns.
For i = 1, ..., n, let yi ∈ Rc be the vector corresponding to the i-th row of U . We cluster the points
(yi )y=1,...,M in Rc with the k-means algorithm into clusters C1 , ..., Cc , where c was determined by
visual inspection of the grouped elements of W after clustering (see Fig.1).
2.3

Personalized multi-task neural network for pain level estimation

As in [3], we use a multi-task neural network (MT-NN) approach with shared layers and task-specific
layers. While in [3] each task corresponded to a different person, here we assign a task to each profile.
The benefits of using profiles as tasks are two-fold: (i) more data is available to tune the models and
avoid over-fitting, and (2) in real-world applications, when a new subject arrives, only a small amount
of data will need to be acquired to assign the subject to a profile, without the need to train a new
layer in the MT-NN. Specifically, in this work our regression MT-NN consisted of one shared hidden
layer, and one task-specific layer. For all units in the MT-NN, we employed the rectified linear unit
(ReLU) activation function: xi = max(0, Mi xi + bi ), where xi represents the input of the i-th layer
of the network (x0 is the input feature vector), and Mi and bi are the weight matrix and bias term.
3

Model
Kächele et al. [19]
Physiology (NC)
Video (NC)
Multimodal (NC)
Multimodal (c=2)

Multimodal (c=3)
Multimodal (c=4)

Cluster 1
Cluster 2
All
Cluster 1
Cluster 2
Cluster 3
All
All

MAE
0.99
1.00
0.92
0.88
0.86
0.77
0.82
0.86
0.73
0.84
0.80∗
0.77∗

RMSE
1.16
1.28
1.23
1.21
1.22
1.16
1.19
1.22
1.11
1.19
1.17∗
1.15∗

ICC(3,1)
N/A
0.28
0.21
0.29
0.22
0.37
0.30
0.19
0.42
0.27
0.32
0.31

Table 2: Model performance in
terms of mean absolute error (MAE),
root-mean-square error (RMSE) and
intra-class correlation (ICC) using
different modalities and clustering
approaches (NC for no clustering
(c=1), and c={2,3,4} for the clustered MT-NNs). Two-tail t-tests
were performed on the "all" conditions against "multimodal (NC)"
and ∗ indicate significance with
p-value ≤ 0.05.

The proposed model was implemented using deep learning frameworks TensorFlow 1.2.1 [22] and
Keras 2.0.6 [23]. To optimize the network parameters, we first trained a joint network (c=1) using the
Adam optimizer [24], with mean absolute error loss. We then fixed the weights of the shared layer,
and used the learned weights of layer i = 2 to initialize the profile-specific layers, which were further
fine-tuned using the data corresponding to the subjects assigned to each profile. To regularize the
network, we applied dropout [25] and employed an early stopping strategy based on a validation set.

3

Results

For each recording, we used the first part, corresponding to the first 48 pain stimuli, as training
set, the second part corresponding to the following 10 stimuli as validation set, and the final part
corresponding to the final 22 stimuli as test set. The training set was used to cluster each of the
85 subjects (2 subjects were excluded due to poor landmark tracking by OpenFace) into c clusters
or profiles. Several meta-information statistics of these clusters were calculated. They are shown
in Table 1 and indicate that our clustering process does not result in any differences in the gender
composition of the clusters, but it discriminates according to the subject’s age.
Once the cluster assignment was completed, we extracted features from overlapping windows as
described in Sec.2.1 We then balanced the training set to have equal amount of P = 0 and P > 0
instances by downsampling the over-represented class (P = 0). As performance measures, we
used the mean absolute error (MAE), root-mean-square error (RMSE), and intra-class correlation
ICC(3,1) [26]. The results are shown in Table 2, and also compared to related work [19], which used a
multimodal early fusion of geometric and appearance based features, SC, ECG, and electromyography
(EMG). Our results indicate that similar performance is achieved with the proposed multi-modal
(unclustered) approach and the prior work. The results also indicate an overall improvement of the
multimodal approach with respect to the single-modality approaches, and also by the models with
the proposed clustered multi-task approach. Furthermore, we note that the profiling approach also
finds clusters in which pain estimation performance is better, whereas other clusters seem to be more
challenging. This can be seen from the scores of cluster 3 in the MT-NN with c = 3, where ICC
reaches 42% and achieves the best MAE and RMSE errors, thus indicating that the MT-NN model
achieves improved estimation performance on this subgroup. However, it drops in performance on
the other two clusters, which shows the need to further investigate the data of those subjects. By
comparing the models with different number c of profiles, on average we obtain similar performance,
with c = 4 performing the best in terms of MAE and RMSE, while c = 3 attains the best ICC score.

4

Conclusions

We proposed a clustered multi-task neural network model for continuous pain intensity estimation
from video and physiological signals. Each task in our model represents a cluster of subjects with
similar pain response profiles. We showed the benefit of our multimodal multi-task model with
respect to (a) single-task (population) models, and (b) single-modality approaches. We conclude that
the choice of cluster profiles would need to be selected based on target application and output metric.
Future work will focus on improving the profiling approach and optimizing the network topology for
exploiting their commonalities.
4

References
[1] Amanda C. de C. Williams and Kenneth D. Craig. Updating the definition of pain. PAIN, 157(11):2420–
2423, 11 2016.
[2] Jarred Younger, Rebecca McCue, and Sean Mackey. Pain outcomes: A brief review of instruments and
techniques. Current Pain and Headache Reports, 13(1):39–43, 2009.
[3] Daniel Lopez-Martinez and Rosalind Picard. Multi-task Neural Networks for Personalized Pain Recognition
from Physiological Signals. In Seventh International Conference on Affective Computing and Intelligent
Interaction Workshops and Demos (ACIIW), San Antonio, TX, 2017.
[4] Roi Treister, Mark Kliger, Galit Zuckerman, Itay Goor Aryeh, and Elon Eisenberg. Differentiating between
heat pain intensities: The combined effect of multiple autonomic parameters. Pain, 153(9):1807–1814,
2012.
[5] Daniel Lopez Martinez, Ognjen Rudovic, and Rosalind Picard. Personalized Automatic Estimation of
Self-Reported Pain Intensity from Facial Expressions. In 2017 IEEE Conference on Computer Vision and
Pattern Recognition Workshops (CVPRW), pages 2318–2327, Hawaii, USA, 7 2017. IEEE.
[6] Dianbo Liu, Fengjiao Peng, Andrew Shea, Ognjen Rudovic, and Rosalind Picard. DeepFaceLIFT:
Interpretable Personalized Models for Automatic Estimation of Self-Reported Pain. IJCAI 2017 Workshop
on Artificial Intelligence in Affective Computing, 2017.
[7] Audrey-Anne Dubé, Marco Duquette, Mathieu Roy, Franco Lepore, Gary Duncan, and Pierre Rainville.
Brain activity associated with the electrodermal reactivity to acute heat pain. NeuroImage, 45(1):169–180,
3 2009.
[8] Kenneth M. Prkachin. Assessing pain by facial expression: Facial expression as nexus. Pain Research and
Management, 14(1):53–58, 2009.
[9] Kenneth M. Prkachin and Patricia E. Solomon. The structure, reliability and validity of pain expression:
Evidence from patients with shoulder pain. Pain, 139(2):267–274, 2008.
[10] M J L Sullivan, D A Tripp, and D Santor. Gender differences in pain and pain behavior: The role of
catastrophizing. Cognitive Therapy and Research, 24(1):121–134, 2000.
[11] Steffen Walter, Sascha Gruss, Hagen Ehleiter, Junwen Tan, Harald C Traue, Stephen Crawcour, Philipp
Werner, Ayoub Al-Hamadi, and Adriano O Andrade. The biovid heat pain database data for the advancement and systematic validation of an automated pain recognition system. In 2013 IEEE International
Conference on Cybernetics (CYBCO), pages 128–131. IEEE, 6 2013.
[12] Mathias Benedek and Christian Kaernbach. Decomposition of skin conductance data by means of
nonnegative deconvolution. Psychophysiology, 47:647–658, 2010.
[13] D.M. Alexander, C. Trengove, P. Johnston, T. Cooper, J.P. August, and E. Gordon. Separating individual
skin conductance responses in a short interstimulus-interval paradigm. Journal of Neuroscience Methods,
146(1):116–123, 7 2005.
[14] W. Engelse and C. Zeelenberg. A single scan algorithm for QRS detection and feature extraction. Computers
in Cardiology, 6:37–42, 1979.
[15] Philipp Werner, Ayoub Al-Hamadi, Robert Niese, Steffen Walter, Sascha Gruss, and Harald C. Traue.
Automatic Pain Recognition from Video and Biomedical Signals. In 2014 22nd International Conference
on Pattern Recognition, pages 4582–4587. IEEE, 8 2014.
[16] Tadas Baltrusaitis, Peter Robinson, and Louis-Philippe Morency. OpenFace: An open source facial
behavior analysis toolkit. In 2016 IEEE Winter Conference on Applications of Computer Vision (WACV),
pages 1–10. IEEE, 3 2016.
[17] Erroll Wood, Tadas Baltruaitis, Xucong Zhang, Yusuke Sugano, Peter Robinson, and Andreas Bulling. Rendering of Eyes for Eye-Shape Registration and Gaze Estimation. In 2015 IEEE International Conference
on Computer Vision (ICCV), volume 2015 Inter, pages 3756–3764. IEEE, 12 2015.
[18] Paul Ekman and Wallace Friesen. Facial action coding system. 2002.
[19] Markus Kächele, Mohammadreza Amirian, Patrick Thiam, Philipp Werner, Steffen Walter, Günther Palm,
and Friedhelm Schwenker. Adaptive confidence learning for the personalization of pain intensity estimation
systems. Evolving Systems, 8(1):71–83, 3 2017.
[20] Jianbo Shi and Jitendra Malik. Normalized cuts and image segmentation. IEEE Transactions on Pattern
Analysis and Machine Intelligence, 22(8):888–905, 2000.
[21] Ulrike von Luxburg. A tutorial on spectral clustering. Statistics and Computing, 17(4):395–416, 12 2007.

5

[22] Martín Abadi, Paul Barham, Jianmin Chen, Zhifeng Chen, Andy Davis, Jeffrey Dean, Matthieu Devin,
Sanjay Ghemawat, Geoffrey Irving, Michael Isard, Manjunath Kudlur, Josh Levenberg, Rajat Monga,
Sherry Moore, Derek G Murray, Benoit Steiner, Paul Tucker, Vijay Vasudevan, Pete Warden, Martin Wicke,
Yuan Yu, and Xiaoqiang Zheng. TensorFlow: A System for Large-Scale Machine Learning. 12th USENIX
Symposium on Operating Systems Design and Implementation (OSDI 16), pages 4265–283, 2016.
[23] François Chollet et al. Keras, 2015.
[24] Diederik P. Kingma and Jimmy Ba. Adam: A Method for Stochastic Optimization. In Proceedings of the
3rd International Conference on Learning Representations, 2015.
[25] Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov. Dropout:
A Simple Way to Prevent Neural Networks from Overfitting. Journal of Machine Learning Research,
15:1929–1958, 2014.
[26] Patrick E. Shrout and Joseph L. Fleiss. Intraclass correlations: Uses in assessing rater reliability. Psychological Bulletin, 86(2):420–428, 1979.

6

