On Time Domain Analysis of Photoplethysmogram Signals for Monitoring Heat Stress

There are a limited number of studies on heat stress dynamics during exercise using the photoplethysmogram (PPG) and its second derivative (APG). However, we investigate the most suitable index from short PPG signal recordings for heat stress assessment. The APG waveform consists of a, b, c and d waves in systole and an e wave in diastole. Our preliminary results indicate that the use of the energy of aa area, derived from PPG signals measured from emergency responders in tropical conditions, is promising in determining the heat stress level using 20-s recordings. After examining 14 time domain features using leave-one-out cross-validation, we found that the aa energy extracted from PPG signals is the most informative feature for classifying heat-stressed subjects, with an overall accuracy of 79%. Moreover, the combination of the aa energy with the traditional heart rate variability index of heat stress (i.e., the square root of the mean of the squares of the successive aa intervals) improved the heat stress detection to an overall accuracy of 83%.

However, this is not a feasible option for heat stress detection due to variations in skin sweating. Therefore, PPG measurement that can be measured in an ambulatory setting represents an attractive option. As Rowell et al. [6] found reductions in cardiac output, central blood volume and stroke volume during heat stress, we expect that the morphology of the PPG waveform will change accordingly. Recent improvements in wearable sensor technology allow for the continuous measurement of PPG for the purpose of measuring heart rate. However, information derived from the PPG signal can go beyond exclusive heart rate measurement to provide additional information on physiological and autonomic arousal and stress conditions.
Our recent investigation to detect systolic peaks in PPG signals measured after exercise showed some challenges due to motion artifacts, sweat and non-stationary effects [7]. Several filters and algorithms have been examined to analyze of the PPG wave contour; however, they still lack accuracy and reproducibility [8]. As a result of these challenges, researchers have started to apply the second derivative to emphasize and easily quantify the delicate changes in the PPG contour [9]. For these reasons, a second derivative will be used in this project to increase accuracy.
The acronym for the second derivative of the PPG is usually SDPPG or APG; however, APG will be used exclusively in this study, according to the recommendations in [10]. As shown in Figure 1b, the APG waveform consists of four systolic waves (a, b, c, and d waves) and one diastolic wave (e wave) [11]. The height of each wave was measured from the baseline, with the values above the baseline being positive and those under it negative. 15   Fingertip photoplethysmogram signal morphology [12]. (a) Fingertip photoplethysmogram; (b) second derivative wave of photoplethysmogram. The photoplethysmogram waveform consists of one systolic wave and one diastolic wave, while the second derivative photoplethysmogram waveform consists of four systolic waves (a, b, c, and d waves) and one diastolic wave (e wave).
The relative heights of these waves (b/a, c/a, d/a and e/a ratios), particularly the b/a ratio, are related to aging and carotid distensibility [13] and are used in calculating the aging index (b − c − d − e)/a [14]. Recently, the detection of a waves in APG signals has been used to calculate heart rate [15,16] and heart rate variability indices [17][18][19]. Despite the application of APG to cardiac variables and that the clinical significance of APG measurement has been well investigated, there is still a lack of studies focusing on heat stress assessment using APG signals. Matsuyama [20] assessed the most suitable index for heat stress assessment in APG signals, determining that the typical indicators of APG analysis, such as c/a, d/a and e/a, are insufficient for heat stress assessment, as they are very difficult to detect correctly: c, d and e waves merge after heat stress [21]. Therefore, our investigation is focused on examining time domain features based on systolic references in heat-stressed PPG and APG signals, especially after the automatic detection of a and b waves became feasible [22].

Data Collection
The heat stress PPG data for this study were collected as part of a National Critical Care and Trauma Response Centre (NCCTRC) project to assess the physiological and perceptual responses of emergency responders to a simulated chemical, biological and radiological (CBR) incident in tropical environmental conditions, to compare the efficacy of various cooling methods. The background of the NCCTRC's thermal research can be found in [23]. Forty healthy, heat-acclimatised emergency responders (30 males and 10 females) with a mean ± SD age of 34.7 ± 6.6 volunteered and provided written informed consents to participate in this study, which was approved by the Human Research Ethics Committee of the Northern Territory Department of Health and Menzies School of Health Research. Participants undertook 30 min of triaging and resuscitating, transporting and decontaminating weighted manikins, while wearing Level 3 personal protective equipment, which comprised a fully-enclosed, impermeable suit, including, boots, gloves, hood, face mask and respirator (SE400i, S.E.A. Group, Warriewood, Australia), followed by 30 min of rest and cooling. This protocol was repeated three times with PPG data collected during each rest period, as shown in Figure 2. The database is available upon request at Charles Darwin University: http://www.cdu.edu.au/ehse.  Here, PPG data were measured by a photoplethysmography-equipped device (Salus APG, Osaka, Japan) at a sampling rate of 367 Hz, with the sensor located at the cuticle of the second digit of the left hand. Measurements were taken for 20 s while participants were undertaking seated rest. An emergency physician annotated the systolic peaks as controls for evaluation. The participants were normotensive (mean systolic blood of 129.3 mmHg, range 110-165 mmHg) and had no known cardiovascular, neurological or respiratory disease. Prior to the experiment, the participants provided information about their physical condition. Physical information, such as height and weight, were also measured for reference and summarized in Table 1. Alcohol consumption and smoking were eliminated during 24 h and 2 h before experimentation, respectively. For signal analysis, MATLAB 2010b (The MathWorks, Inc., Natick, MA, USA) was used. An Omron HEM-907 was used for blood pressure measurement.

Methodology
In this study, our goal is to find out the most informative feature that can be used for heat stress assessment. Furthermore, we investigate whether the PPG or APG waveform potentially improves the heat stress analysis.

APG Signal
To produce an informative APG signal, two steps are required: filtering the PPG and differentiating the filtered PPG. We applied a zero-phase second-order Butterworth filter, with a 0.5-7 Hz bandpass based on the optimization process in [22]. Then, the second derivative was applied to the filtered PPG to analyze the APG waveform. Equations (1) and (2) represent a non-causal filter; the three-point center derivative was created with a delay of only two samples, as follows: where T is the sampling interval and equals the reciprocal of the sampling frequency and n is the number of data points. Figure   (i) One heart beat of the PPG and its APG measured from a subject at rest; (ii) one heart beat of the PPG and its APG measured from a subject after simulated heat stress induction.

Feature Extraction
To our knowledge, this is the first study to investigate the correlation of PPG features with the effect of heat stress. Therefore, there was no feature tested before to capture the heat stress impact on PPG signals. To start this line of research, it was logical to examine existing features used in the literature to diagnose different diseases, such as the b/a index, the amplitude of the a wave and the amplitude of the b wave. Furthermore, we tested new features to determine the optimal feature for heat stress detection, such as the energy of the aa area, the energy of the ab area, the energy of the ba area and the slope of the ab segment. In total, we tested 14 time domain features, seven features extracted from the PPG signals and seven features extracted from the APG signals, which can be defined in a compact manner as follows: (3) where S refers to the either PPG or APG signal, N refers to the total number of beats in the processed signal and A and B arrays contain the annotated a waves and b waves, respectively, in each APG signal. Colon notation is used to specify sequences, for example S(x : y) is a vector of PPG or APG data points from x to y.
As the heart rate may significantly increases due to heat stress, we calculated beat-to-beat intervals from the APG signals using a waves [16]. The annotated a waves are used to calculate the duration of each consecutive aa interval as follows: where A contains the annotated a waves in each APG signal, n is the number of a waves and aa contains the a-a intervals.

Statistical Analysis
We calculated seven features from each PPG and APG signal recording. As we had 40 subjects, each feature set contains 80 values. Each feature set consisted of 40 values calculated from subjects measured before the simulated heat stress induction and 40 values calculated from subjects measured after the simulated heat stress induction. As we have a small sample size, we compared the values within each feature set by applying the paired Mann-Whitney test (p < 0.05 was considered significant). Because we considered all of these features simultaneously, it is likely that a few p-values are small merely due to stochastic fluctuations rather than due to systematic differences between subjects measured before and after the simulated heat stress induction. As a consequence, the p values need to be appropriately corrected. One may try to control the probability that a false positive occurs by applying a Bonferroni post-correction [24]. As we are dealing with many different simultaneous tests (42 tests in total), it is more natural to try to control the false discovery (false positive) rate. Therefore, we used the Holm-Bonferroni method because it controls the false positive rate and is a simple test that is uniformly more powerful than the Bonferroni correction [25].
As the main objective of this research is to find the optimal time domain feature for assessing simulated heat stress induction, we tested four classifiers: Mahalanobis distance, linear discriminant analysis (LDA), quadratic discriminant analysis (QDA) and the linear support vector machine (SVM). The classification rates were calculated using leave-one-out cross-validation. Two statistical measures were used for the output of each classifier: sensitivity (SE), which was calculated using the formula TP/(TP + FN), and positive predictivity (PP), which was calculated using the formula TP/(TP + FP); whereas TP is the number of true positives (subjects measured after simulated heat stress induction detected in heat-stressed subjects); FN is the number of false negatives (subjects measured after simulated heat stress induction have not been detected as heat-stressed subjects); and FP is the number of false positives (subjects measured before simulated heat stress induction detected as heat-stressed subjects). To compare the performance of the features given the small dataset in each classifier, we applied the F 1 score, as recommended in [26], which is defined as 2 × (SE × PP)/(SE + PP). The overall accuracy (OA) is the average of F 1 scores for all classifiers.

Results and Discussion
An ingestible telemetric temperature sensor (CorTemp 100, HTI Technologies, Florida, MI, USA) was used to measure core body temperature by transmitting a signal proportional to the temperature of the gastrointestinal tract to a handheld receiver. Once located in the gastrointestinal tract, the thermosensitive pill is a valid index of core temperature when referenced to rectal or esophageal temperature [27,28].
Participants ingested a core temperature pill with fluid prior to breakfast that preceded the commencement of data collection by approximately 4 h. This time frame was adopted to allow the pill to empty from the stomach and enter the small intestine while minimizing the risk of the pill being emptied from the body prior to the completion of the study [27]. To control for the possibility of local cooling of the telemetry pill via ingestion of fluids and/or crushed ice, subjects were excluded from the study if gastrointestinal temperature was less than 35.5 • C and/or decreased by 2 • C in any 5-min period during the rest phase. Figure 4 demonstrates a significant difference between at-rest measurement and those after exercises based on the BCT. We tried to collect the PPG signals immediately after exercise. However, some participants may have cooled down while queuing for measurement. Therefore, the obtained p-value from comparing at-rest measurement and the simulated heat stress simulations did not linearly decrease. For example, the returned p-value from the first BCT test (at rest vs. after exercise 1) was 1.6 × 10 −13 , while the second test (at rest vs. after exercise 2) was 3.3 × 10 −4 and for the last test (at rest vs. after exercise 3) was 0.0036. Note, this step was carried out to ensure heat stress impact.
It is expected that the aa interval will be shorter as the heart rate increases during exercise. Figure 5 illustrates the significant difference between the four measurements based on the average aa intervals per participant using the paired Mann-Whitney test. It is clear that the aa interval shows an overall significant difference (between at-rest measurement and those after exercise) compared to the BCT. The returned p-value from the first aa test (at rest vs. after exercise 1) was 9.1 × 10 −8 , while the second test (at rest vs. after exercise 2) was 3.6 × 10 −8 and for the last test (at rest vs. after exercise 3) was 3.9 × 10 −8 . However, the aa interval will be used in calculating the RMSSD. Therefore, we excluded the aa interval as a time domain feature from our investigation.  Figure 4. Boxplot of body core temperature (BCT) for subjects measured before and after exercise. The p-value was calculated using the paired Mann-Whitney test.
The separability performance of the discussed 14 time domain features for predicting heat stress is shown in Table 2. It is clear that the energy of the ba area of the PPG signal provides the lowest overall p-value (4.1 × 10 −9 ), while the energy of the aa is the second lowest overall p-value (4.8 × 10 −9 ). The lowest and second lowest features are very close in terms of separability. As the main goal of this study is to find the optimal feature with significant separability power, we examined the classification performance of each feature using four different classifiers, Mahalanobis, LDA, QDA and SVM, as shown in Table 3. Table 2. Feature comparison in differentiating between before and after exercise measurements. Mean and standard deviation values of 14 features extracted from the PPG and APG waveforms. The p-value of discriminating between before exercise (BE) and after three exercises (E1, E2 and E3) is calculated using the Mann-Whitney test. The average of the three p-values (p) is given in the last column. Uncorrected p-values from the Mann-Whitney test, where * and ** indicate p < 0.05 and p < 0.005, respectively; † indicates p values that remain significant with a p < 0.05 after post-correction (Bonferroni-Holm, α < 0.05). 3.6×10 5 ± 4.7×10 5 6.4×10 5 ± 6.6×10 5 1.
After applying the leave-one-out cross-validation, the energy of aa scored the highest overall accuracy (OA = 79.40%), while the energy of ba scored the second highest overall accuracy (OA = 77.88%). Now, after using four different classifiers in all heat stress simulation stages, we have clear evidence that the energy of aa is the optimal feature for detecting heat stress. Figure 6 shows the LDA classification performance of the optimal feature for detecting heat stress after three exercises. The optimal feature (energy of aa) detected heat stress after the first exercise with a SE of 82.5% and PP of 75%, as shown in Figure 6a. The PP was slightly improved in classifying after exercise 2 (cf. Figure 6b); however, the detection of heat stress was considerably improved after exercise 3 with an SE of 87.5% and a PP of 79.55%, as shown in Figure 6c. Interestingly, the overall accuracy of detecting heat stress increased as the exercise progressed from one simulated heat stress induction stage to the next. This new outcome demonstrates that the energy of the aa area can be used to assess heat stress without the need for measuring the heart rate, which can be used as an indicator of heat stress tolerance.     (c) at rest vs. after exercise 3. Here, RMSSD stands for square root of the mean of the squares of the successive aa (heartbeats) intervals; SE stands for sensitivity; and PP stands for positive predictivity. The plus signs refer to subjects measured at rest, while the diamond signs refer to subjects measured after the simulated heat stress induction.
To ensure that the energy of the aa area can be used to replace other measures of heat stress, it is important to compare its performance with another heat stress index, such as BCT and RMSSD. The problem is we could not measure BCT immediately after exercise, and a few subjects had to wait, which cooled them down a bit (cf. Figure 4). Therefore, we had to compare the optimal feature with the RMSSD calculated from the PPG signals. Figure 7 shows the performance of the RMSSD for detecting all simulated heat stress inductions. The RMSSD detected heat stress after the first exercise with an SE of 77.5% and a PP of 58.49%, as shown in Figure 7a. The detection of heat stress was considerably improved after exercise 2 and exercise 3 with an SE of 87.5% and a PP of 63.64%, as shown in Figure 7b,c. It is clear that the optimal feature (the energy of aa) outperformed the RMSSD heat stress index. Note that the RMSSD showed no classification improvement in detecting heat stress after exercise 2 and Exercise 3: the sensitivity and positive predictivity remained the same (cf. Figure 7b,c).
To maximize the heat stress detection performance, we combined the optimal feature (energy of aa) and the traditional HRV heat stress index (RMSSD). Figure 8 shows the performance of the combined features (energy of aa and RMSSD) for detecting all simulated heat stress inductions. The RMSSD detected heat stress after the first exercise with an SE of 75% and a PP of 75%, as shown in Figure 8a. The detection of heat stress was considerably improved after exercise 2, with an SE of 85% and a PP of 79.1%, as shown in Figure 8b. The progression of detecting heat stress continued, and the accuracy improved after exercise 3 with an SE of 92.5% and a PP of 86.1%, as shown in Figure 8c. It is clear that the use of combined features noticeably improved the heat stress detection performance with each sequential heat stress simulation induction stage.

Limitations of the Study and Future Work
Exploring these findings across different genders, ages and health condition groups will improve the generalization across the population. To our knowledge, there is no available PPG database measured in tropical conditions or after heat stress that would allow a more thorough assessment and comparison of the tested algorithms. In future studies, it may be advisable to have multiple PPG systems to collect the signal immediately after exercise. In the present study, data were collected immediately after exercise. However, some participants may have cooled down while queuing for measurement.

Conclusions
The findings of this preliminary study indicate that heat stress can be assessed using the blood volume of PPG signals. This new outcome shows that combining the energy of the aa area with the RMSSD of the PPG signals can be used to assess heat stress and can be used as an indicator of heat stress tolerance. The results of this study indicate that PPG can be a potential modality for heat stress analysis and identification of individuals at risk. Our preliminary study demonstrates indicative results and now motivates the need for a larger study that validates the energy of the aa area from the PPG against traditional heat stress tolerance indices. The developed classifier can then be embedded in wearable PPG sensors to ultimately detect the stress level in real time.