Very Short-Term Photoplethysmography-Based Heart Rate Variability for Continuous Autoregulation Assessment

Featured Application

A signal processing combination of instantaneous pulse rate variability and time shift multiscale entropy for autoregulation assessment was proposed.

Abstract

Background: Heart rate variability (HRV) has been widely applied for disease diagnosis. However, the 5 min signal length for HRV analysis is needed. Method: A signal processing procedure for very short-term photoplethysmography (PPG) signal for fever detection and autoregulation assessment was proposed. The Time-Shift Multiscale Entropy Analysis (TSME) was applied to instantaneous pulse rate time series (iPR) and normalized by the cumulative distribution function (CDF) of all scales to calculate novel indices. A total of 33 subjects were recruited for the study. Fifteen participants whose body temperatures were higher than 37.9 °C were served as the fever group. Others were served as the non-fever group. The total 15 s PPG signal with 200 sampling rates was used for iPR calculation. Result: The CDF value of entropy on the scale k = 19 (CDF(E(k = 19))) of iPR had the lowest p-value calculated by the Weltch t-test between two groups (p < 0.001). The Spearman correlation r between CDF(E(k = 19)) and body temperature is −0.757, 0.287, and −0.830 in all subjects, the non-fever group and the Fever group, respectively. The area under the curve, calculated from the receiver operating characteristic of CDF(E(k = 19)) of iPR is 0.915. Conclusion: The entropy of iPR is useful for detecting fever. Moreover, a short-term PPG signal is suitable to develop real-time applications, and multiscale entropy provides different scales of information for daily healthcare.

1. Introduction

Heart rate variability (HRV) has been applied for health assessment for a long time. Many studies have shown that HRV is related not only to cardiovascular diseases but also to autonomic nervous system (ANS) activity [1,2]. To measure HRV, the beat-to-beat interval (RRi) of the electrocardiography (ECG) signal is usually calculated. In time-domain analysis, indices such as Standard deviation of normal -to-normal intervals (SDNN), Proportion of NN50 divided by the total number of NNs (pNN50) calculated from RRi are useful for the diagnosis [3]. In frequency domain analysis, each frequency band in the power spectrum of RRi has been proved related to the ANS activity [2]. Furthermore, the irregularity of RRi in abnormal heartbeats is different from healthy heartbeats [4,5,6]. Instead of using ECG for HRV calculation, photoplethysmography (PPG) signal can also be used for HRV analysis, called pulse rate variability (PRV). PPG was usually worn on the user’s finger. PPG optically measures the change in blood volume in the microvascular tissue. The pulse-to-pulse interval (PPi) is obtained for PRV analysis. Although pulse arrival time (PAT) exists between the ECG and PPG signals, previous studies still showed that PRV can be a surrogate for HRV [7,8,9].

Both RRi and PPi have low time resolution, which limits the investigation of HRV and PRV [1,10]. Information about every two heartbeats is usually ignored. Therefore, instantaneous pulse rate variability (iPRV) was proposed [11]. iPRV calculates the instantaneous frequency of the main intrinsic mode function (IMF), decomposed by empirical mode decomposition (EMD) [12]. According to Bedrosian’s theorem, the source data should be no amplitude modulation for instantaneous frequency calculation. Therefore, the instantaneous frequency is calculated using the normalized Hilbert transform (NHT) or the normalized direct quadrature (NDQ) [13]. For comparison of RRi or PPi, the inverse of instantaneous frequency, called instantaneous pulse rate (iPR), is usually applied for further analysis. Time-domain and frequency-domain analyses of iPRV have been studied [11,14,15]. These studies showed that iPRV had a similar result to conventional indices compared with HRV or PRV in time- and frequency-domain analysis. Moreover, because of the higher time resolution of time series iPR, more detailed information on the variability of heartbeats can be indicated. However, the non-linear dynamical analysis of iPRV is not performed yet.

For nonlinear dynamical analysis, entropy is common for complexity measurement of time series. Shannon, who proposed Shannon entropy, was the first person who applied the concept of entropy to information theory [16]. However, Shannon entropy is only related to the probability of the appearance of the element, and it cannot measure the complexity of randomness. In the field of biomedical signal analysis, sample entropy (SampEn) is the most popular one for nonlinear dynamical analysis [17]. SampEn is based on conditional entropy, which compares the difference of short patterns in the signal. For HRV analysis, RRi was applied SampEn for complexity measurement. Many studies show that the RRi in patients with cardiovascular disease were different from the RRi in healthy people [5,6,18]. In addition, the different time scales of the signal are also concerned for more detailed information indication [5,6]. In this study, the multiscale entropy method called Time-Shift Multiscale Entropy analysis (TSME) [19] was applied in the iPR to indicate the difference between the person with fever and the person without fever. Body temperature and heart activity are usually regulated by ANS and peripheral response [20]. Fever is one of the early signs of infection and needs to be systematically monitored in the intensive care unit (ICU) [21]. Furthermore, fever is also a common sign of many infectious diseases and non-infectious diseases, and it is also an indicator of disease progression in many hospitalized patients. Body temperature measurement using traditional thermometers is reliable and low cost. However, it requires caregivers to perform the measurement and it will be time-consuming to measure frequently. Therefore, it is hard for clinicians to continuously monitor patients’ body temperature and get notified if a patient has a fever in time using the traditional technique. PPG had been widely used in the hospital setting to measure oxygen saturation and estimate blood pressure. Obtaining temperature information from the PPG can assist the clinicians to have a better understanding of the patients’ condition. In addition, PPG devices have been installed on several models of smart wearing devices on the market nowadays. Being able to monitor fever using PPG in an out-of-hospital setting may allow healthy people to get aware of their illness earlier, such as COVID-19 or heat stroke during exercise, and have a chance to seek medical assistance in time.

The structure of the paper is organized as follows: Section 2 describes the experiment, subject information, and algorithms of each signal processing and analysis method. Section 3 gives the result of TSME on iPR and PPi. The discussion of the results and further comparison are described in Section 4. Finally, Section 5 gives a conclusion.

2. Materials and Methods

2.1. Subjects and Experiment

This experiment was approved by the Research Ethics Committee for the Protection of Human Subjects of the National Chiao Tung University (NCTU-REC-103-061; date of approval: 2 June 2015). Informed consent was obtained from all subjects and their parents before the experiment. The experiment was carried out at the Yo Yo Clinic, Kaohsiung, Taiwan. A total of 33 people aged 7 to 18 years were recruited for this study. The subject’s core body temperature was measured by the doctor. 15 subjects whose body temperature was higher than 37.9 °C and who were fever diagnosed by the doctor served as the fever group, others served as the non-fever group. Subjects were asked to do supine position quietly for ten minutes. At the same time, the PPG signal was collected with 200 sampling rates by using the PPG device (Nonin 8500, Nonin Medical Inc., Plymouth, MN, USA). The previous study showed that iPRV can be performed using short-term PPG [22], a total of 15 s stable signal in the middle of the total PPG signal was used in this study. It can also decrease the computation time cost of the advanced EMD method and the SampEn algorithm.

2.2. PRV Procedure

The PRV procedure is similar to HRV but uses the PPG signal instead of the ECG signal. The first step of PRV is to find the pulse peak of the PPG signal. Then calculate the time difference between a pulse peak and the next one. The PPi can be obtained after calculating the interval between each two pulse peaks. However, in this study, only 15 s PPG was applied. The time scale of PPi was not enough for further analysis. Therefore, the PPi was interpolated into 200 sampling rates, which was the same as the PPG signal to be PPi series for further analysis.

2.3. iPRV Procedure

The iPRV applied the main decomposition method EMD of the Hilbert-Huang Transform (HHT) to calculate the instantaneous pulse rate of PPG [11]. IMF was decomposed by EMD from PPG. The IMF should satisfy two conditions: (1) the number of local minimum and local maximum must be equal to the number of zero-crossings or differ at most by one, and (2) the value of average envelop must be approximately equal to zero at any point. However, the mode mixing problem, which means that different frequency band components are decomposed into the same IMF, exists in IMFs if the signal is by EMD. Therefore, the complete ensemble EMD (CEEMD) was proposed to solve the problem [23]. Before introducing CEEMD, EMD should be described first. The first step of EMD is to calculate the upper envelope and the lower envelope by using extremum points of the raw signal. Next, compute the average envelope using the upper envelope and lower envelope. Third, subtract the average envelope from the raw signal as the output signal. Finally, check if the output signal is IMF or not. If not, let the output signal be a raw signal to redo the above steps. If the output signal is IMF, output the IMF, and using the raw signal subtract the IMF to redo the above steps until the residual is mono-component. For CEEMD, the raw signal is added and subtracted by the same white noise before performing EMD. Therefore, there are two types of raw signal: raw signal with adding white noise and raw signal with subtracting white noise. Decompose these two raw signals into IMFs. To decrease the influence of white noise, repeat the above steps with different white noise many times (50 times in this study). Finally, average all corresponding output IMFs to ensemble IMFs. According to Nuttall’s theorem, before computing the instantaneous frequency by the Hilbert transform or direct quadrature method, the oscillation of the source data should be a sinusoidal signal. That is why we decomposed PPG by EMD into IMFs. Additionally, according to Bedrosian’s theorem, the source data should be without amplitude modulation for instantaneous frequency calculation by Hilbert transform or the direct quadrature method. Therefore, in this study, the NDQ was applied to calculate the instantaneous frequency [13]. After the modulation of the amplitude was removed, the source data can be considered as a cosine function. Therefore, we can use the direct quadrature method to calculate the instantaneous phase using the arctangent. The instantaneous frequency is obtained from the derivative of the instantaneous phase. For comparison with the RRi or PPi series, the inverse of instantaneous frequency called the instantaneous pulse rate series (iPR) was computed for further processing.

2.4. Time-Shift Multiscale Entropy (TSME)

The TSME applied the main idea of Higuchi’s fractal dimension (HFD) to compute the different time scales of the source data [19,24]. Let X be the source data.

$$X=x_1,x_2,x_3,\dots ,x_N$$

(1)

The HFD reconstructs the source data into different time scales as follows:

$$y_k^{\beta }=\left(x_{\beta },x_{\beta +k},x_{\beta +2k},\dots ,x_{\beta +{\displaystyle \lfloor }\frac{N-\beta }{k}{\displaystyle \rfloor }k}\right),\beta \in Z^{+},0<\beta <\frac{N-\beta }{k}$$

(2)

where N is data length, $y_k^{\beta }$ is kth time scale segment. $x_{\beta }$ is the βth point of source data. The TSME of the kth time scale is calculated:

$$TSM{E}_k=\frac{1}{k}{\displaystyle \sum }_{\beta =1}^{k}SampEn(y_k^{\beta })$$

(3)

In this study, the TSME were m = 2, r = 0.15* standard deviation of the iPR or PPi series, k = 1 to k = 20. The E(k = i) denotes the ith TSME value. Figure 1 showed the demonstration of signal processing. Figure 2 showed the example of the distribution of E(k = 1~20) of a subject. Further indices according to the distribution were calculated. The TSME of all k scales could be normalized by cumulative distribution function (CDF). The CDF value of TSME in scale i is denoted as (CDF(E(k = i))).

2.5. Other Indices Calculation

In addition to the above entropy indices, the other entropy-related indices and conventional iPRV and PRV indices were calculated. Mean(E(k = 1~20)) and std(E(k = 1~20)) calculated the average value and standard deviation of E(k = 1~20). Max(E(k = 1~20)) and Min(E(k = 1~20)) was the maximum value and the minimum value of E(k = 1~20). Area(E(k = 1~20)) calculated the distribution area of distribution of E(k = 1~20).

The conventional time domain indices were the mean and standard deviation of time series (iPR or PPi). For the frequency domain analysis, the frequency bands of low frequency (LF), high frequency (HF), and very high frequency (VHF) were 0.04 Hz to 0.15 Hz, 0.15 Hz to 0.4 Hz, and 0.4 Hz to 0.9 Hz, respectively. The normalized power of each band was calculated as follows.

$$\mathrm{nLF}, \mathrm{or} \mathrm{nHF}, \mathrm{or} \mathrm{nVHF}=\frac{\left(\mathrm{LF}, \mathrm{or} \mathrm{HF}, \mathrm{or} \mathrm{VHF}\right)}{\mathrm{TP}}$$

(4)

where TP is the total power. However, VHFs were not used in conventional PRV analysis. Only the iPRV analysis used VHF. Therefore, the TP can be calculated including VHF or not. Thus, we use superscript a to present the normalized power calculated excluding VHF and superscript b to present the normalized power calculated including VHF.

3. Results

The basic characteristic of the subjects is shown in

Table 1. Basic information of subjects.

	Non-Fever Group	Fever Group	Total
# of subjects	18	15	33
# of boy	11	10	21
age	10.78 ± 3.15	10.4 ± 2.38	10.60 ± 2.79
Body temperature (°C)	36.28 ± 0.30 *	38.4 ± 0.47	37.25 ± 1.14

* p < 0.05 compared to the Fever group.

.

The result related to iPR is shown in Figure 3. TSME on the scale k = 19 CDF(E(k = 19)) had the lowest p-value calculated using the Weltch t-test between two groups (p < 0.001).

The Spearman correlation r between CDF(E(k = 19)) and body temperature is shown in Figure 4. The r values are −0.757, 0.287, and −0.830 in all subjects, the non-fever group and the fever group, respectively.

The AUC, calculated from the ROC, of each entropy variables are shown in

Table 2. The AUC of variables related to the entropy value of iPR.

Variable	AUC	Sensitivity	Specificity
Mean(E(k = 1~20))	0.659	0.667	0.722
std(E(k = 1~20))	0.681	0.733	0.611
Max(E(k = 1~20))	0.567	0.467	0.778
Min(E(k = 1~20))	0.657	0.667	0.778
Area(E(k = 1~20))	0.663	0.667	0.722
E(k = 4)	0.685	0.733	0.778
CDF(E(k = 19))	0.915	0.933	0.833

The AUC of CDF(E(k = 19)) had the highest value of 0.915, the sensitivity is 0.933, and the specificity is 0.833. Figure 5 showed the AUC distribution of E(k) and CDF(E(k)) calculated by using iPR.

For comparison with conventional PRV time series,

Table 3. The AUC of the variables related to the entropy value of PPi.

Variable	AUC	Sensitivity	Specificity
Mean(E(k = 1~20))	0.678	0.600	0.722
std(E(k = 1~20))	0.689	0.800	0.667
Max(E(k = 1~20))	0.678	0.800	0.611
Min(E(k = 1~20))	0.678	0.800	0.556
Area(E(k = 1~20))	0.678	0.600	0.722
E(k = 4)	0.689	0.800	0.611
CDF(E(k = 19))	0.745	0.733	0.722

showed the AUC of the above variables calculated by using the PPi series. The AUC of all variables related to the entropy value of the PPi series was lower than the CDF(E(k = 19)) of the iPR.

4. Discussion

4.1. Body Temperature and iPRV

Since humans are warm-blooded animals, we need to maintain our body temperature in a fixed range. However, while our bodies are infected by germs, to increase the power of immune body temperature, the corresponding organs should be regulated by ANS [20,25]. Increase heart rate by the sympathetic nervous system (SNS) and vasoconstriction, for example. Some studies showed the relationship between thermoregulation and HRV [26,27]. Their result showed that the very low frequency band power (VLF) had a corresponding change, while the core body temperature or the environment temperature changed. However, they did not discuss the situation of fever. Moreover, to indicate a meaningful VLF, the signal length should be longer [1]. Our previous studies still showed that conventional nLF and nHF in patients with fever were different from those without fever [28]. Furthermore, we also examined whether the conventional frequency band indices in iPRV are similar to the indices in conventional HRV in our previous work in the time and frequency domain [14,15]. nLF and nHF calculated by iPRV actually showed the difference in ANS activity between patients with fever and people without fever. For non-linear dynamic analysis, the higher resolution time series can present much detail about the complexity of the time series. The CDF(E(k = 19)) of the iPR series also had a higher AUC than the CDF(E(k)) of PPi, whatever k. Therefore, the following term CDF(E(k = 19)) is denoted as the one of iPR series.

Table 4. The AUC of time and frequency domain indices of iPRV.

Variable	AUC	Sensitivity	Specificity
Mean(iPR)	0.914	0.867	0.889
std(iPR)	0.730	0.733	0.722
LF	0.7	0.667	0.722
HF	0.826	0.733	0.889
VHF	0.763	0.733	0.778
LF/HF	0.620	0.667	0.722
nLF a	0.620	0.667	0.722
nHF a	0.620	0.667	0.722
nLF b	0.630	0.533	0.944
nHF b	0.672	0.8	0.556
nVHF b	0.507	0.733	0.444
Pulse rate	0.944	0.933	0.889

^a calculated total power including VHF; ^b calculated total power excluding VHF.

and

Table 5. The AUC of time and frequency domain indices of PRV.

Variable	AUC	Sensitivity	Specificity
Mean(PPi)	0.944	1	0.889
std(PPi)	0.792	0.8	0.833
LF	0.733	0.8	0.667
HF	0.837	0.867	0.889
LF/HF	0.702	0.733	0.722
nLF	0.702	0.733	0.722
nHF	0.702	0.733	0.722

show the AUC calculated by time- and frequency domains of iPRV and PRV. Only pulse rate and mean (PPi) had a higher AUC than CDF(E(k = 19)). In addition to Mean(PPi), Mean(iPR) also had a high AUC, but it is because they are directly correlated with pulse rate. Although pulse rate differences are found between patients with fever and non-fever (p < 0.001), only pulse rate to find patients with fever is not enough. Influencing pulse rate is so easy. For example, the pulse rate of someone after exercise also possibly had a high pulse rate which is similar to the pulse rate of patients with fever. The pulse rate distribution is shown in Figure 6a. The pulse rate in the fever group is significantly higher than in the non-Fever group. The difference between pulse rate and CDF(E(k = 19)) is that CDF(E(k = 19)) had a high correlation compared to body temperature in the ever group, but the correlation was low between pulse rate and body temperature in the fever group (Figure 4c and Figure 6d). For patients suffering from fever, autoregulation of not only their ANS but also peripheral response. The CDF(E(k = 19)) may correspond to this reflection of autoregulation or the severity of the disease. CDF(E(k = 19)) had a higher absolute correlation |r| compared to pulse rate in each group and total subjects (Total: 0.757 vs. 0.726; Non-fever group: 0.287 vs. 0.05; Fever group: 0.830 vs. 0.267). It showed that CDF (E (k = 19)) is not directly affected by pulse rate. Body temperature, ANS activity, or peripheral response may be the possible factors that affect the CDF(E(k = 19)), especially in the Fever group. Applying logistic regression using pulse rate and CDF(E(k = 19)) can increase AUC. Using these two indices to detect patients with fever may be more reliable (

Table 6. Logistic model of pulse rate and CDF(E(19)).

Logistic Model	AUC	Sensitivity	Specificity
−1.28 + 0.02 Pulse rate-1.57CDF(E(k = 19))	0.952	0.933	0.889

and Figure 7).

4.2. Apply Short-Term PPG for CEEMD and Entropy Analysis

Figure 8 showed the time cost of SampEn, CEEMD, and the total time cost. The total time cost is dominated mainly by SampEn. When the sample length is larger than 3000 points, the total time cost is around 2 to 4 s. In addition, five minutes of heartbeats is needed for conventional HRV analysis [1]. These limitations restrict the development of real-time HRV application. Therefore, in this study, only 15 s were recorded for signal processing. However, the output of CEEMD had a boundary effect problem, the first 2.5 s signal and the last 2.5 s signal were eliminated. Total 10 s signal used for iPRV and TSME analysis. The time cost is only about 2 s. Because iPRV breaks the time resolution problem in HRV, the resolution of short-term time series is still enough for time and frequency-domain analysis [11,22]. For non-linear dynamic analysis, the SampEn is sensitive to sample length [17]. Some research suggests that the sample length N should be larger than 200 points for more consistent findings [29]. However, the effect of multiscale entropy analysis still needs further examination. In addition, the basic algorithm of SampEn, the cost of time is O(n²), was applied in this study [17]. Some algorithms can reduce the time cost of SampEn calculation [30]. Moreover, the parameter settings of SampEn and TSME are important. In this work, we follow the same setting as Pham, who proposed the TSME [19]. However, it is unknown whether it is the best setting for this study. In future studies, the testing of a different combination of parameter settings is necessary. Research proposed the quadratic sample entropy (QSE) method to dynamically vary the parameter r to reduce the influence of the parameter [31]. Applying the QSE may be a solution for further study.

4.3. Fever Detection and Application

To measure body temperature, the most common and easy way is using a thermometer. However, it is hard to measure body temperature continuously with a thermometer. Another way to measure body temperature is to use a thermographic camera. The thermography camera is usually used for measuring the body temperature measurement at the same time. However, the price of a thermographic camera is not cheap. Therefore, in this study, we proposed a new method for detecting fever by using a very short-term PPG signal. Although our method cannot measure body temperature directly, detecting fever is more important. Especially in the ICU, it is necessary to monitor patients who have a fever. Fever is also a common sign of many infectious diseases and non-infectious diseases. Using a thermometer to measure body temperature is a cost of manpower and cannot be measured continuously. Using a thermographic camera is too expensive. However, PPG has been widely used in ICU for pulse oximetry and breathing measurement. Moreover, measuring PPG through the camera of a mobile is also available. In addition, iPRV can also provide not only ANS but also peripheral circulation information using short-term PPG [11,14,22,28]. For homecare or other application, PPG is a cheap and non-invasive measurement, and can be applied to many situations. For example, PPG can be integrated with the watch or the camera of a cell phone can measure the PPG signal as well. Furthermore, the non-contact PPG is also developing [32,33,34,35]. The PPG has the potential to be the health monitor in many situations.

4.4. Limitations

Some limitations remain in this study. The first is a small population of subjects. The second is that only 7 to 18 years are recruited for this study. Although the response to autoregulation in 7 to 18 years is stronger, whether the age influence the result of this study is unknown. Therefore, to recruiting more subjects and researching the influence of age is further work. The third is we only recruited the fever as a high body temperature group. Fever is not the only way to increase body temperature. For example, someone after exercise also may have a high body temperature. The difference between the fever group and the after-exercise group depends on whether the set point changes or not. In other words, if a person without fever has a high body temperature, the hypothalamus and ANS will tend to decrease body temperature. For example, a decrease in heart rate or vasodilation. Current research has examined the difference between the fever group and the low body temperature group. However, it is currently unknown whether research can distinguish the reasons for high body temperature is unknown. Recruiting other groups with high body temperature due to other reasons to compare the difference among each other is also important for further research.

Additionally, the PRV result is only for reference in this study. Actually, for conventional PRV or HRV analysis, the 5 min signal is suggested [1]. Nevertheless, for comparison with iPRV, this study only applied a 15 s PPG signal for PRV analysis which may lead to the result being meaningless or not stable enough. However, the result of the iPR series or iPRV would not be influenced by this limitation.

5. Conclusions

In this study, we propose a new signal processing method for continuous fever detection by using a very short-term PPG signal. The result showed that the multiscale entropy of the iPR series is useful for fever detection. Moreover, the variables calculated by iPRV and TSME had the potential to be an indicator of autoregulation. The short-term PPG signal is suitable to develop real-time application and multiscale entropy provide different scales of information for daily healthcare.