Nonlinear Heart Rate Variability features for real-life stress detection. Case study: students under stress due to university examination

1Stress has been investigated as a risk factor for cardiovascular disease and for reduced human performances, which in some situation, such as dangerous works or driving a car, may results in negative consequences. Stress influences the balance of Autonomous Nervous System (ANS).

2HRV is a non-invasive measure reflecting the variation over time of the period between consecutive heartbeats (RR intervals) and has been proved to be a reliable marker of ANS activity.

3For this reason, several studies investigated cardiovascular reaction induced by stress using Heart Rate Variability (HRV) focussing on acute, laboratory stressors: cognitive (e.g., mental arithmetic), psychomotor (e.g., mirror tracing) challenges and physical stressors. Moreover, as standard laboratory stressors do not always engage subjects' affective response, real life stressors (e.g. precompetitive anxiety or social interaction stressors such as public speaking tasks) are often applied to provide a more appropriate social context in which negative emotions might be elicited. Some studies investigated HRV variations in the case of university exams as it is a real-life stressor. These studies included only linear HRV measurement, except for the study by Anishchenko which considered nonlinear measures such as Approximate Entropy. In the current study, we investigated how the most common nonlinear HRV measures vary in subject under stress due to university examination. Furthermore, we proposed a classifier for automatic detection of stress based on nonlinear HRV features.

4We performed a prospective analysis, examining 5-minute HRV extracted from ECG records of volunteer students in two different conditions: the first record was performed during an on-going verbal examination (stress session); the second one was performed after holidays (control session).

5The data were acquired from 42 students of the School of Biomedical Engineering of the University Federico II, who volunteered to take part in the study. This study was performed in compliance with the Human Study Committee regulations of the University of Naples "Federico II". After obtaining written consent, a 3-lead electrocardiogram (ECG) was recorded on 2 different days: the first recording was performed during an ongoing university verbal examination (stress session), while the second one was taken in controlled resting condition (rest session) after a holiday period, far away from stress induced from study routines.

6There are many factors that may influence the HRV, such as circadian rhythm, body position, activity level prior to recording, medication, verbalization and breathing condition. For that reason, we took special precautions to maintain similar condition, such performing both recordings at similar time of day and in a sitting body position after an adaptation time of at least 15 minutes. Furthermore, we asked about consumptions of drugs, and none of the students declared consumption of drugs. Finally, we induced participants to speak also in the control session.

7We performed a short-term 5-minute HRV analysis according to International Guidelines. The RR interval time series were extracted from ECG records using an automatic QRS detector, WQRS available in the PhysioNet's library, based on nonlinearly scaled ECG curve length feature. Two scientists independently reviewed and corrected the QRS detection and manually labelled the normal beats obtaining the so called series of normal to normal (NN) beat intervals. QRS review and correction was performed using PhysioNet's WAVE. The fraction of total RR intervals labelled as normal-to-normal (NN) intervals was computed as NN/RR ratio. This ratio has been used as a measure of data reliability, with the purpose to exclude records with a ratio less than a 90% threshold. None of the records were excluded as NN/RR is higher than 90%.

8Nonlinear properties of HRV were analyzed by the following methods: Poincaré Plot, Approximate Entropy, Correlation Dimension, Detrended Fluctuation Analysis, and Recurrence Plot. We focussed on these methods as they were implemented in a software freely distributed and widely used for research activities.

9The Poincaré Plot (PP) is a common graphical representation of the correlation between successive RR intervals, for instance the plot of RRj+1 versus RRj. A widely used approach to analyze the Poincaré plot of RR series consists in fitting an ellipse oriented according to the line-of-identity and computing the standard deviation of the points perpendicular to and along the line-of-identity referred as SD1 and SD2, respectively.

10Approximate entropy measures the complexity or irregularity of the RR series. The algorithm for the computation of Approximate Entropy was briefly described here.

11Given a series of N RR intervals, such as RR1, RR2,..., RRN, a series of vector of length m X1,X2,..., XN-m+1 is constructed from the RR intervals as follows: Xi,=[RRi, RRi+1... RRi+m-1].

12The distance d[Xi,Xj] between vectors Xi and Xj is defined as the maximum absolute difference between their respective scalar components. For each vector Xi, the relative number of vectors Xj for which d[Xi,Xj]≤r, is computed where r is referred as a tolerance value (see equation 1).

13Then, the following index Φm(r) is computing by taking natural logarithm of each and averaging them over i.

14Finally, the approximate entropy is computed as:

15In this study, we computed the ApEn with m = 2 and with three different value of the threshold r:

16r = 0.2*SDNN (standard deviation of the NN series);

17r = rmax that is, the value of r in the interval (0.1 * SDNN, 0.9 * SDNN) which maximizes the ApEn;

18r = rchon, that is the value computed according to the following formula proposed by Chon:

19where N denotes the length of the NN sequence, and SDDS and SDNN, respectively, are the measure of the short-term and long-term variability of the RR sequence. Formally, SDDS is the Standard Deviation of the Difference Sequence of the series RR, that is, [RRi+1 - RRi, RRi+2 - RRi+1,..., RRN - RRN-1], and SDNN is the Standard Deviation of the series NN.

20The value of the parameters r and m were chosen according to the recommendation for slow dynamic time series, such as heart rate variability, (m = 2 and r = 0.2*SDNN) and to the findings of recent studies which suggested choosing the value of r which maximizes the entropy (r = rmax) and proposed a formula for automatic selection of the value r (r = rchon).

21Further in the paper, we will indicate the Approximate Entropy computed with the different values of r with the following notation En(0.2), En(rmax) and En(rchon).

22The correlation dimension D2 is another methods to measure the complexity used for the HRV time series.

23As for Approximate Entropy, the series of Xi is constructed and is computed as in formula 2, but the distance function, in this case, is defined as follows:

24where Xi(k) and Xj(k) refer to the k-th element of the series Xi and Xj, respectively.

25Then, the following index Cm(r) is computed by averaging over i.

26The correlation dimension D2 is defined as the following limit value:

27In practice this limit value is approximated by the slope of the regression curve (logr,logCm(r)). In the current study a value of m = 10 was adopted.

28Detrended Fluctuation Analysis measures the correlation within the signal and consists into the steps described here.

29The average of the RR interval series is calculated on all the N samples. The alternate component of RR interval series, which is defined as RR minus its average value, is integrated:

30The integrated series is divided into non-overlapping segments of equal length n. A least square line is fitted within each segment, representing the local trends with a broken line. This broken line is referred as yn(k), where n denotes the length of each segment.

31The integrated time series is detrended as following: y(k)-yn(k). The root-mean-square fluctuation of the detrended time series is computed according to the following formula:

32The steps from 2 to 4 are repeated for n from 4 to 64.

33Representing the function F(n) in a log-log diagram, two parameters are defined: short-term fluctuations (α1) as the slope of the regression line relating log(F(n)) to log(n) with n within 4-16; long-term fluctuations (α2) as the slope of the regression line relating log(F(n)) to log(n) with n within 16-64.

34Recurrence Plot (RP) is another approach performed for measurement of the complexity of the time-series. RP was designed according to the following steps.

35Vectors X i= (RRi, RRi+τ,..., RRi+(m-1) τ), with i = 1,..., K, with K=[N-(m-1)* τ)], where m is the embedding dimension and τ is the embedding lag, are defined.

36A symmetrical K-dimensional square matrix M1 is calculated computing the Euclidean distances of each vector Xi from all the others.

37After choosing a threshold value r, a symmetric K-dimensional square matrix M2 is calculated as the matrix whose elements M2(i,j) are defined as:

38The RP is the representation of the matrix M2 as a black (for ones) and white (for zeros) image.

39In this study, according to, the following values of the parameters introduced above were chosen:, with SDNN defined as the standard deviation of the NN series.

40In the RP, lines are defined as series of diagonally adjacent black points with no white space. The length l of a line is the number of points which the line consists of.

41The following measures of RP were computed: recurrence rate (REC) defined in equation 11; maximal length of lines (lmax); mean length of lines (lmean); the determinism (DET) defined in equation 12; the Shannon Entropy (ShEn) defined in equation 13.

42The HRV analysis was performed using Kubios for all the measures except the Approximate Entropy ones which were computed using in-house software in Matlab as their computation is not available in Kubios. All the computed measures are summarized in Table 1.

43We calculated mean, standard deviation, median and 25th and 75th percentiles to describe distribution of HRV features during stress and rest conditions. Moreover, we calculated mean, standard deviation, median and 25th and 75th percentiles of the individual differences between stress session and rest session, and we used the Wilcoxon signed rank test to investigate the statistical significance of features' variation within each subject. The statistical analysis was performed by in-house software developed in MATLAB version R2009b (The MathWorks Inc., Natick, MA).

44We adopted Linear Discriminant Analysis (LDA) as classification method. LDA aims to find linear combinations of the input features that can provide an adequate separation between two classes, in the current study, stress and rest session. LDA uses an empirical approach to define linear decision plans in the feature space. The discriminant functions used by LDA are built up as a linear combination of the variables that seek to maximize the differences between the classes. Further details about LDA can be found in Krzanowski.

45In order to evaluate the classifier, we computed the common measures for binary classification performance measurement using the formulae reported in Table 2, considering positive to the test those records classified as under stress. Total classification accuracy represents the ability of the classifier to discriminate between the two sessions, sensitivity refers to the ability to identify records in the stress session and specificity refers to the ability to identify records in the rest session.

46To estimate the performance measures we adopted a 10-fold cross-validation scheme. This technique consists in developing 10 classifiers as following: (1) dividing randomly the dataset into 10 subsamples; (2) excluding a subsample (testing subset) in turn; (3) developing a classifier with the remaining 9 subsamples (training subset); (4) testing each classifier with the excluded subsample (which is used as an independent testing dataset), computing the performance measures using the formulae in Table 2. The 10-fold cross-validation estimates of the performance measure are computed as the averages over the 10 classifiers. We divided the dataset in 10 folds by subject and not by record in order to obtain a person-independent testing.

47It would be possible to use all the 13 selected HRV features reported in Table 1 for the LDA, however this may decrease the performance of the classifier, particularly because of curse of dimensionality. Therefore, we tried to find the subset of features which could discriminate the two classes with the highest total classification accuracy: we adopted the so-called exhaustive search method, investigating all the possible variations with repetition of k out of N features (with k from 1 to N). Since the number of features N is 13, we investigated 213 = 8192 subsets of features, training and testing the same number of classifier, as discussed in the previous subsection.

48For all the single features and for the best subset of features, that is, which achieved the highest total classification accuracy, the discrimination function was computed against all the dataset in order to provide classification rules.

49All the analysis was performed by in-house software developed in MATLAB version R2009b (The MathWorks Inc., Natick, MA).

50Table 3 shows the descriptive statistics of nonlinear short-term HRV measures in the enrolled subject during the rest and stress sessions. Table 4 presents how nonlinear short-term HRV measures vary in the subjects in rest or under stress due to an ongoing university examination. The features SD2, D2, En(0.2), En(rchon), α1, lmax were significantly reduced during university examination as compared with rest session, while lmean, REC and ShEn increased significantly during stress.

51Table 5 shows the performance of the classifier based on single nonlinear HRV measures. The highest 10-fold cross-validation estimate of the total classification accuracy was achieved by SD2 and by D2. Applying the linear discriminant analysis against the whole dataset, we obtained the following rules: for instance, referring to SD2, if SD2 is lower than 0.0646 ms, the record is classified as under stress, otherwise as in rest condition. Moreover, also the classifiers based on REC, En(rchon) and α1 obtained a satisfying accuracy rate (71%).

52The classifier achieving the highest accuracy is based on the subset of features SD1, SD2 and En(0.2), obtaining a total classification accuracy rate of 90%. All the performance measures are reported in the last row of Table 5. The classification rule can be express as follows:

53The record is classified as stress if:

54Furthermore, the classification rule could be represented as in Figure 1: in the 3D space of the features SD1, SD2 and En(0.2), the points under the plan defined by equation (12), are classified as STRESS; those above as REST.

55Among the classifier based on couple of features for comparison with other studies it is interesting to report the performance of the classifier based on SD1 and SD2 which achieved a total classification accuracy, a sensitivity and a specificity rate of 82%, 79% and 86%.

56In this study, we compared within-subject variations of short-term nonlinear HRV measures in healthy subjects during condition of mental stress due to an on-going university examination.

57Almost all the features measuring complexity of the time series statistically decreased during the stress session, like D2, En(0.2), En(rchon), which have been widely used complexity measures for HRV.

58Almost all the features measuring complexity of the time series statistically decreased during the stress session. These findings confirms the results obtained by Anishchenko, which showed that Approximate Entropy decreased significantly during stress condition due to university examination. Among the approximate entropy measures considered in this study, the one based on threshold value rchon achieved the highest total classification accuracy. These results support previous findings regarding rchon capability to detect different physiological conditions.

59Furthermore, our finding of decreased complexity measures, in particular D2, are in line with studies about the relationship between Heart Rate complexity and acute physical stress or short-term mental stress.

60The decreased value of complexity measures reflects a change towards more stable and periodic behaviour of the heart rate under stress which may be associated with stronger regularity, decoupling of multimodal integrated networks and deactivation of control-loops within the cardiovascular system. As interpreted by Schubert, this reduction in heart rate complexity during a high stress condition may reflect a lower adaptability and fitness of the cardiac pacemaker.

61The results of the classification for automatic detection of high-stress reinforce the findings of the statistical analysis: the D2 and the En(rchon) enables detecting the stress condition with a total classification accuracy rate higher than 70%. Furthermore, also the SD2, which is a measure of long-term variability, and α1, which provided information about short-term fluctuations, achieved comparable performances.

62The combination of features achieving the best results consists of the two parameters of PP (SD1 and SD2) and a measure of complexity (En(0.2)) and enables detecting the stress condition with a total classification accuracy, a sensitivity and a specificity rate of 90%, 86% and 95%, respectively. The SD1 was chosen in the best combination of features, although the classifier based only on SD1 achieved the lowest performance among the one based on single features, because it provided information different from the other features, particularly SD2 and En(0.2).

63We underlined that, even if not shown in best combination, the classifier based on the two parameters of the PP (SD1 and SD2) achieved a total classification accuracy higher than 80%, confirming the usefulness of PP as a valid marker for mental stress.

64The performance achieved by the selected subset of nonlinear features is higher than that achieved by selected linear feature on the same data-set reported in our unpublished observations. Furthermore, comparing with the study of Kim, who adopted a logistic regression on linear HRV features for distinguishing high stressed subject from low stressed ones, achieving a total classification accuracy of 70%, the performance of the current study are better. These comparisons confirms the usefulness of nonlinear HRV features for automatic classification.

65Controlled breathing was not asked in order not to affect student performance during the university exam. However, the effect of breathing pattern on HRV is a debated question. Some studies showed that different breathing conditions may have an impact on the reproducibility of HRV. In contrast, other studies found that such factors did not have a significant impact on HRV reliability and their findings seem to suggest that HRV is reliable and consistent over time, whether or not respiration is controlled.

66In the current study, we focussed only on a few nonlinear methods, those which were implemented in Kubios, a free software for HRV analysis. Although this choice is a limit of the current study, it could be useful in order to increase the reproducibility of the experiment by other investigators.

67As regards the classification methods, LDA succeeded partially in separating the two classes, providing an intelligible model. The intelligibility of features and classification rule is strongly appreciated in medical domain data-mining. However, the adoption of a linear classifier may represent another limit of the current study, which did not enable us to consider nonlinear structures in classification. In future work we will use nonlinear methods such as Artificial Neural Network (ANN) and Support Vector Machine (SVM) with adequate kernel, in order to achieve a possible improvement in the performance measurement. However, we underlined that the computational cost of the LDA is lower than the ANN or SVM, saving time in the operation.

68Finally, the results of this paper could extend the use of portable sensing devices, usually adopted in cardiac applications, to stress detection.

69In conclusion, the results of the current study suggest that nonlinear HRV analysis using short term ECG recording could be effective in automatically detecting real-life stress condition, such as a university examination. The proposed classifier based on the Poincaré Plot measures and on the Approximate Entropy enables detecting the condition of stress due to university examination with a total classification accuracy, a sensitivity and a specificity rate of 90%, 86, and 95%, respectively.

70Further research on a large sample size and on different stressful conditions will help to further elucidate the findings of this study and effectiveness of HRV analyses for differentiation between low and high stress condition.

71Abbreviation

72ACC: Total classification accuracy; ANS: Autonomous Nervous System; ECG: Electrocardiogram; HRV: Heart Rate Variability; LDA: Linear Discriminant Analysis; NPV: Negative Predictive Value; PP: Poincaré Plot; PPV: Positive Predictive Value; RP: Recurrence Plot; SEN: Sensitivity; SPE: Specificity.

73The authors declare that they have no competing interests.

74PM and LP conceived the study, evaluated the data, performed data analyses. PM recruited subjects, managed data acquisition, and drafted the manuscript. MB and LP supervised the study. All authors read, reviewed and approved the final manuscript.