Beat-to-beat continuous blood pressure estimation with optimal feature set of PPG and ECG signals using deep recurrent neural networks

Data and feature extraction

For evaluating the proposed algorithm, we utilized raw data sourced from the MIMIC II dataset^[12]. This dataset comprises a comprehensive collection of photoplethysmogram (PPG) and electrocardiogram (ECG) signals, along with corresponding blood pressure (BP) measurements, obtained from thousands of patients in intensive care units (ICUs). The inclusion of such diverse and extensive data from real-world ICU settings enhances the robustness and generalizability of our algorithm's evaluation. This dataset contains PPG and ECG signals with corresponding BP from thousands of patients in intensive care units (ICU). However, ECG and PPG signals are instable, which can be affected by physiological activities, electromagnetic environment, and other factors. Signal preprocessing is necessary in order to extract features from ECG and PPG accurately. In this case, Kachuee et al. have made available preprocessed and clean data derived from the MIMIC II dataset^[13]. Noises from all signals were removed using an averaging filter^[13]. The dataset, consisting of 943 subjects and containing a total of 12,000 records, is available for access in the UCI machine learning repository. All PPG, ECG, and BP signals are sampled at 125 Hz and only the records longer than 8 minutes are chosen to develop our proposed RNN model. In our work, we used 660 subjects and each subject has one record. DBP and SBP were obtained from the minimum and maximum values in each heartbeat of each subject. Figure 2 illustrates the distribution of diastolic blood pressure (DBP) and systolic blood pressure (SBP). The dataset of 660 subjects was divided into four groups based on the office blood pressure threshold, as detailed in Table 1. Within our dataset, the standard deviation of SBP is 19.47 mmHg, while the standard deviation of DBP is 9.64 mmHg.

Figure 2. SBP and DBP ranges in the dataset.

In order to confine the values of all samples within the range of (0, 1) for each record, it is essential to perform data normalization^[14]. The equation indicated as (1) presents the expression used for normalizing the data. In this equation, the minimum and maximum values of the signals, denoted as x_minand x_max respectively, are obtained from the signals x(t) within a specific window. The resulting normalized signals are represented by $$\tilde{x}(t)$$.

In our work, six types of features with a total number of 28 were extracted from every ECG and PPG beat, as shown in Table 2. In order to extract these features, some fiducial points of each ECG and PPG beat are located, such as R-peaks of ECG beats, pulse foot, systolic peaks, dicrotic notches, and diastolic peaks of PPG beats. Figure 3 displays the schematic diagram illustrating the fiducial points and various features. The R-peaks of the ECG beats are extracted using the hierarchical clustering method^[15] and all fiducial points in PPG beats are extracted based on the wavelet transform method^[16]. In terms of phase space features of PPG signals, phase space reconstruction (PSR) and box-counting techniques are used^[14]. Figure 4 provides an example illustrating how the phase-space diagram is divided into N × N pixels, where N represents an integer value.

Figure 3. Some features of ECG and PPG signals.

Figure 4. PSR features of a PPG beat.

As an illustration, the phase-space diagram of a single PPG beat is plotted using a time delay of 20 milliseconds and N is set to 20. The phase-space diagram consists of pixels, with certain pixels being considered as black boxes (B_all) if they are traversed by one or more trajectories, and the remaining pixels are designated as white boxes (W_all). The count of black boxes (B_all) among all the pixels provides an indication of the spread of these trajectories. Subsequently, the phase-space diagram is divided into nine areas. The complexity or fractal dimension of a specific area in the PSR diagram is defined as d_i= B_i/B_all, i = 1, 2, ..., 9. The index i corresponds to a specific area within the diagram. Here, we only measured the degree of complexity in eight areas because there are no trajectories in area 9.

Feature selection

The efficiency of the model will be affected when a large number of features are used, because some features might be related to others, leading to redundant computational costs. Therefore, feature selection is used to construct the model by selecting the most relevant features among all features. It can reduce the dimensionality of input data of the model^[17]. The objective in this case is to minimize the computational cost of the model while enhancing its performance through the incorporation of the most relevant features. In our work, five feature selection methods will be used, which are named ReliefF feature selection (ReliefF), correlation-based feature selection (CFS), mutual information method (MI), neighborhood component analysis (NCA), and relativity analysis (RA).

Vote system

To address the limitations of feature selection algorithms, a voting system is implemented, as illustrated in Figure 5. First, for each specific algorithm, the importance of each feature was given, which is called the weight of features and features were sorted based on the weights at the same time. Then, for each algorithm, the top N features were assigned one vote and the final voting score sorted to select features. Here, N represents the number of features that are selected to build the model. Finally, features with high correlation will be removed based on the Pearson correlation coefficient.

Figure 5. Flowchart of the proposed vote system.

Deep recurrent neural network

The proposed deep recurrent neural network (RNN) model, as depicted in Figure 6, provides an overview of its architecture. It consists of a Bi-LSTM at the bottom layer, which is followed by a stack of five LSTM layers with residual connections. Furthermore, at the top two layers, there is a Bi-LSTM layer and a LSTM layer.

Figure 6. The proposed deep RNN architecture. Each circle represents a LSTM cell. Each red dashed box is a Bi-LSTM layer, which consists of forward (green) and backward (blue) LSTM. The blue dashed box represents the LSTM layer with residual connections.

In general, RNNs with deeper architecture tend to outperform shallow RNNs by a significant margin^[25,26]. According to the experimental findings reported in^[26], it was observed that the utilization of multiple layers of RNN led to a notable enhancement in expressive capacity. Nevertheless, employing a fully deep network can result in an escalation of computational complexity. Moreover, deep RNNs can be challenging to train due to issues such as gradient exploding and vanishing. To overcome these problems while preserving the benefits of recurrent connections, specific strategies are employed. An independently RNN (Ind-RNN) model was proposed^[26]. This model effectively addresses the problems of gradient vanishing and exploding by regulating the gradient backpropagation through time.

An independently RNN (Ind-RNN) model was proposed^[26]. It can regulate the gradient backpropagation through time to effectively address the gradient vanishing and exploding problems. The degradation problem was effectively addressed in^[27] through the introduction of a deep residual learning framework. The experimental results indicated that their proposed deep residual nets achieved significant accuracy improvements as a result of the substantial increase in depth. Inspired by the deep residual learning framework in^[27], we introduce a residual connection block designed to enhance the propagation of gradients within deep RNNs, specifically LSTM-based models. This approach facilitates more efficient gradient propagation in deep RNN architectures. An example of residual connections is given in Figure 7. The residual connections are implemented between a^th LSTM layer and layer (a = 1, 2, 3, …, n). In this case, $$ x_{t}^{a}$$ and $$h_{t}^{a}$$ are the input and hidden state of the a^th LSTM layer at t^th moment, respectively. $$ \mathcal{H}$$ represents the LSTM function. The input of the LSTM layer at t^th moment is then calculated by adding the input and hidden state of the a^th LSTM layer ( $$ x_{t}^{a}$$ and $$ h_{t}^{a}$$). The integration of residual connections within the LSTM block can be accomplished by utilizing the following equations:

Figure 7. An example of LSTM with residual connections.

Finally, the proposed deep RNN model can be developed by sequentially stacking multiple LSTM blocks of this type. Once the hidden state of the top layer (a = n) of LSTM is calculated, the final output y_t at t^th moment can be computed by equation (19):

In our work, the depth of the teacher model was set as 8. The experimental results demonstrated that when the depth of the deep RNN model exceeds 8, it can lead to the issue of overfitting. The reason for this phenomenon might be the limited training samples of VA classification. In addition, the sequence length of each training sample is initially set to 32, but it can be flexibly adjusted to higher as needed when employing a deeper RNN model. The hidden state size in the deep RNN model is 128 at each layer.

Feature selection

The proposed feature voting system was introduced in Section “METHODOLOGY” in order to select the optimal feature sets for BP estimation. The Pearson correlation coefficient was calculated based on 28 features, setting ρ_x,y < -0.8 or ρ_x,y > 0.8 as strong correlation. Table 3 shows the features with high correlation among 28 features. Finally, 11 different feature sets with the number of features from 7 to 17 are determined in our work, which are shown in Table 4.

Model performance

This section presents the outcomes of our proposed deep RNN model tested on 660 subjects based on different feature sets using 10-fold validation. All features are derived from ECG and PPG signals and each feature set is determined using a vote system of five feature selection methods.

To facilitate comparisons between our proposed model and those found in the literature, we evaluate the performance of blood pressure (BP) estimation using the following equations, which measure estimation accuracy through root mean square error (RMSE) and mean absolute error (MAE).

where y_i is the ground truth of BP and $$\hat{y}_{i}$$ is the BP estimated by the proposed model. Table 5 displays the performance of the proposed model based on different feature sets. MAE and RMSE of SBP, DBP, and MAP are the estimation accuracy metrics. For BP estimating, the best result of the overall performance is achieved using a 14-feature set. The values of RMSE for predicted SBP and DBP are 3.223 and 1.781, respectively, while the MAE values for SBP and DBP estimation are 2.514 and 1.383, respectively. The statistical hypothesis testing is applied on all feature sets via a feature voting system, and a significance level of 0.05 was used. The results are shown in Table 5, where * indicates that the 14-feature set is significantly better than its competitor. Figure 8 shows the performance of DBP and SBP estimation using a 14-feature set. In general, both DBP and SBP are performed in close agreement. We also compared the results with and without containing PSR features. The results show that for both DBP, SBP and MBP, both RMSE and MAE for 14 features with PSR are smaller than those of 14 features without PSR.

Figure 8. Comparison of the ground truth and prediction of a representative subject.

The Bland-Altman graph (BA graph) is widely used in biomedical research to compare the agreement between two measurements^[28]. It shows the agreement between the prediction BP and the reference BP by using the mean difference BP and the average BP^[4]. In Figure 9, the Bland-Altman (BA) graph displays the average value of the reference BP and estimated BP on the x-axis, while the y-axis represents the difference between them. The bias (Mean) and the limits of agreement (Mean ± 1.96 × SD) of our work are shown as dashed line. The mean difference of DBP and SBP are -0.13 and -0.59, respectively. It is shown that for both DBP and SBP, the majority of data points fall within the limits of agreement, signifying a close agreement between the estimated BP and the reference BP. Additionally, Figure 10 illustrates correlation graphs for DBP and SBP. The x-axis represents the ground truth BP, while the y-axis represents the predicted BP. The red dashed line indicates when the ground truth and predicted BP are equal. The green line signifies the linear regression relationship between ground truth BP and predicted BP. On the test data, the correlation coefficients are 0.634 for DBP and 0.888 for SBP.

Figure 9. Bland-Altman plots of the DBP (left) and SBP (right) predictions by the proposed RNN model.

Figure 10. Correlation plot for DBP (left) and SBP (right) targets.