Discrimination of bilateral finger photoplethysmogram responses to reactive hyperemia in diabetic and healthy subjects using a differential vascular model framework

The endothelium is a single layer of cells covering the internal surface of blood vessels and other anatomical structures such as heart valves (Verma and Anderson 2002). By releasing chemical substances such as nitric oxide (NO), the endothelial layer controls blood pressure via changing the vessel diameter and keeps blood coagulation balanced (Sejda et al 2005). Endothelial dysfunction is commonly referred to the disability of this layer to dilate vessels (Anderson et al 1995). It has been established that there is a close relationship between coronary artery disease and endothelial dysfunction even before clinically detectable manifestations of stenosis (Sejda et al 2002). Cardiovascular risk factors such as hypertension, high cholesterol and diabetes cause eventually endothelial dysfunction (Bonetti et al 2003, Iiyama et al 1996, Pellegrino et al 2005), which may in turn be used as a warning sign of early arterial stenosis (Cox et al 1989). Therefore, detecting of the onset of vessel impairments can assist the diagnosis of atherosclerosis at an early stage, making treatment more efficient.

Endothelial function (EF) assessment can be done by using Doppler for blood flow measurement in the coronary artery in response to infusion of vasodilators, which may be pharmacological substances (e.g. nitroglycerin, isosorbate dinitrate, bradykinin or acetylcholine). Since these methods require injection, they are unsuitable for studies involving asymptomatic subjects and multiple measurements for the same individual. It was the introduction of a pioneering technique (Celermajer et al 1992) which offered the advantage of being non-invasive and therefore applicable to EF assessment among asymptomatic subjects. The principle is to expose the endothelium to a physical stimuli (shear stress) using the reactive blood flow to ischemia induced by a pressure cuff (Verma and Anderson 2002). In this method, the brachial artery (BA) flow-mediated dilation (FMD) is evaluated using high-resolution ultrasound imaging during reactive hyperemia, a state of transient increase in organ blood flow that occurs following a brief period of ischemia such as arterial occlusion (Klabunde 2012). FMD is calculated as the percentage of the maximum BA diameter change as compared to the baseline (pre-occlusion) diameter (Corretti et al 2002).

Since the pioneering work of Celermajer, it has been established that cardiovascular events occur more often in patients with impaired FMD of the BA (Neunteufl et al 2000, Gokce et al 2002, 2003), which is positively correlated to coronary artery disease (Takase et al 1998, Zhang et al 2000). Being a systemic disorder (Bonetti et al 2003), impairments in large arteries affect small resistance arteries as well (Iiyama et al 1996, Pellegrino et al 2005, Teragawa et al 2005).

Although non-invasive, the major drawback of the ultrasound-based EF assessment (Celermajer et al 1992) is the requirement for a skilled operator to operate the ultrasonic probe in order to get the appropriate imaging plane and accurately measure the vessel diameter (Corretti et al 2002). Furthermore, ultrasonography can be considered only reliable in assessing the EF in vessels of reasonable size (diameter above 2.5 mm) as smaller vessels are difficult to image (Patel and Celermajer 2006).

Other techniques have been introduced to non-invasively assess the EF by the pulse wave analysis (PWA) of the arterial pressure waveforms (Wilkinson et al 2002), or the digital volume pulse acquired from a photoplethysmograph (Millasseau et al 2000, 2002). In these studies, PWA was combined with the administration of vasodilators such as nitroglycerine. More recently, EF was assessed by comparing the pulse amplitude of the photoplethysmogram (PPG) before occlusion and after cuff release accompanied by ultrasonic FMD (Zahedi et al 2008, Selvaraj et al 2009). Although this method does not require a skilled operator, only a modest correlation (about 55%) was found with ultrasound FMD.

In an attempt to investigate the effect of diabetes on the PPG (Shariati et al 2008), a differential model of vascular system (Zahedi and Alauddin 2004) was utilized. Bilateral PPG signals during rest condition and without any stimuli were recorded, and the results showed a sensitivity of 71.7% and specificity of 70.2%.

More generally, bilateral photoplethysmography has been investigated in multiple studies. Multi-site peripheral PPG (ears, thumbs and big toes) measurement only on healthy subjects showed that the PPG from the right and left sides are highly correlated (more than 98%) at each segmental level (Allen and Murray 2000). Peripheral blood pulsations were also recorded from the fingers and toes of both arms and legs with a multi-channel PPG device in control subjects and patients diagnosed with arterial stenosis in a leg. The results showed a high correlation (Pearson's coefficient r = 0.93) between the bilateral differences in the local blood pressure and in the corresponding pulse wave transit time (PWTT) delay, as well as between the PWTT delay and the ankle-brachial index (r = −0.96) (Erts et al 2005). Finally, PPG signals were recorded bilaterally at the big toe sites for diagnosis and grading of lower limb peripheral vascular occlusive disease using a color relation analysis classifier that resulted in high accuracy and 100% sensitivity (Lin 2011).

In this paper, based on a parametric differential approach (Zahedi and Alauddin 2004), the PPG has been modeled with the objective of discriminating between healthy and diabetic subjects. The model was applied to data recorded from both arms, whereas one arm undergoes cuff occlusion and the other serves as control. The rest of the paper is organized as follows: first, the model, setup and database are described. Then the processing tools are reviewed before results are presented. At the end, the discussion highlights the advantages and limitations of this technique.

2.1. Model

The general structure of this model is shown in figure 1, based on Zahedi and Alauddin (2004). In this linear model, a Gaussian white noise (GWN(n), n as the time index) is the input to the heart, modeled by its transfer function (H_Heart). The output of the heart is a pressure waveform at the aorta, P_Aorta(n). This pressure wave propagates along the vascular tree in upper limbs and is therefore affected by their transfer functions (H_RA,FMD and H_LA). In this study, lower limbs are not considered.

Zoom In Zoom Out Reset image size

Three pressure pulses can be modeled through this model from the same GWN input: P_Aorta, P_LA (LA: left arm) and P_RA,FMD (RA: right arm, under shear stress). Briefly, the PPG-FMD procedure consists of 3 min PPG data recording before occlusion, 4 min during occlusion and 5 min after cuff release (in this work, 3 min of post-release PPG was investigated). Cuff occlusion results in ischemia in the blood-deprived arm, which releases metabolites and reduces the overall hydraulic resistance (Corretti et al 2002). Once the cuff is released, the sudden inflow of blood creates a shear stress which stimulates the endothelial layer (Takase et al 1998, Teragawa et al 2005). After cuff release, the total transfer function of the right arm (H_RA,FMD) is modeled as a combination of endothelial control (H_EC) and controlled transfer function of right-hand arteries (H_RAC). The separation of H_RA,FMD into two parts allows us to emphasize on how does the reactive hyperemia stimulus force the endothelial layer to respond. It should be noted that in this research, H_EC and H_RAC are not treated individually; they are shown with the only purpose of a better understanding of how the FMD test affects the vascular function. Their lumped effect is considered as H_RA,FMD.

Three pressure points indicated in figure 1 denoted by P_Aorta(f), P_LA(f) and P_RA(f) respectively in the frequency domain can be computed by

$$\begin{equation} P_{\rm Aorta} (f) = H_{\rm Heart} (f){\rm GWN}(f) \end{equation} \tag{ 1 }$$

$$\begin{equation} P_{\rm LA} (f) = H_{\rm LA} (f)P_{\rm Aorta} (f)\end{equation} \tag{ 2 }$$

$$\begin{equation} P_{\rm RA} (f) = H_{{\mathop{\rm RA}\nolimits} ,{\rm FMD}} (f)P_{\rm Aorta} (f),\end{equation} \tag{ 3 }$$

where H_Heart, H_LA and H_RA,FMD are the transfer functions of the heart, left arm and right arm respectively (H_RA,FMD: after cuff release). The PPG arises from blood pressure changes at the fingertip. Therefore, the relationship between the pressure waveform and the PPG can be defined by another transfer function, H_P→PPG, and it is assumed that this transformation is the same at both left and right fingers (Reisner et al 2008). Omitting the frequency variable (f) for simplicity and with PPG_LA, PPG_RA,FMD denoting the Fourier transform of PPG_LA(t) and PPG_RA,FMD(t),

$$\begin{equation} {\rm PPG}_{\rm LA} = H_{\rm LA} P_{\rm Aorta} H_{P \to {\rm PPG}} \end{equation} \tag{ 4 }$$

$$\begin{equation} {\rm PPG}_{\rm RA} = H_{\rm RA,FMD} P_{\rm Aorta} H_{P \to {\rm PPG}} .\end{equation} \tag{ 5 }$$

To eliminate the effect of the heart as the main source, a differential approach is applied (Zahedi and Alauddin 2004). The left-arm PPG is considered as the input and the right-arm PPG as the output:

$$\begin{equation} \frac{{\rm PPG}_{\rm RA,FMD}}{{{\rm PPG}_{\rm LA} }} = \frac{{H_{\rm RA,FMD} }}{{H_{\rm LA} }} = H_{{\rm LA} \to {\rm RA,FMD}} . \end{equation} \tag{ 6 }$$

It should be noted that the resulting transfer function (6) between left and right arms would not be physiological. Finally, figure 2 shows the proposed framework in order to illustrate the transition state from the rest PPG signal to the state where the endothelial layer is stimulated using shear stress.

Zoom In Zoom Out Reset image size

2.2. Subjects

2.3. System setup

Two independent, custom-made PPG systems based on commercial modules (Dolphin Medical, Inc.) were utilized to acquire data from the right and left index fingers (for simplicity these two systems have been represented by a single box in figure 3). PPG signals were acquired at a sampling rate of 275 Hz and 16-bit resolution. In order to synchronize the two PPG systems, an external module (Zahedi and Alauddin 2004) was used. The infrared channel of the PPG signal was selected due to its higher signal-to-noise (SNR) power. All data were saved in ASCII form on a personal computer for further analysis.

Zoom In Zoom Out Reset image size

2.4. PPG FMD procedure

2.5. Pre-processing

2.6. System identification

By means of the framework illustrated in figures 1 and 2, system identification techniques (Ljung 1999) were used to extract the model's parameters. The input of the system to be identified consists of the PPG immediately after cuff release from the left index finger (control hand = without endothelial stimuli). The output of this system is the simultaneously PPG recorded from the right index finger (with endothelial stimuli). A fourth-order output-error model structure was found to give satisfactory results, as higher orders resulted in less than 3% improvement on fitness defined by

$$\begin{equation} {\rm Fitness} = 100 \times \left( {\frac{{1 - {\rm norm}(y_m - y_r )}}{{{\rm norm}(y_r - {\rm mean}(y_r ))}}} \right),\end{equation} \tag{ 7 }$$

where y_m is the model output and y_r is the recorded output. The model input–output relationship is given by

$$\begin{equation} y(t) = \frac{{B(q)}}{{A(q)}}u(t) + e(t),\end{equation} \tag{ 8 }$$

where u(t), y(t) and e(t) are the input, output and Gaussian white noise respectively. B(q) and A(q) are the time shift operators defined by

$$\begin{eqnarray} &&B(q) = b_1 + b_2 q^{ - 1} + b_3 q^{ - 2} + b_4 q^{ - 3} \nonumber \\ &&A(q) = 1 + a_1 q^{ - 1} + a_2 q^{ - 2} + a_3 q^{ - 3} + a_4 q^{ - 4} , \end{eqnarray} \tag{ 9 }$$

where a_i and b_i (i = 1,2,3,4) are the parameters of the system. As the response of the endothelial layer to a cuff occlusion stimuli occurs mostly within the first 3 min after releasing the cuff (Zahedi et al 2008, Selvaraj et al 2009), 3 min of immediate post-release left and right PPG were used to identify the model parameters. A total of eight coefficients [b₁, b₂, b₃, b₄, a₁, a₂, a₃, a₄] are thus available as the model representation for each subject. A sample of modeled and measured output is shown in figure 4.

Zoom In Zoom Out Reset image size

2.7. Dimension reduction using singular value decomposition

Singular value decomposition (SVD) has been extensively used for dimension reduction in pattern recognition and information retrieval applications (Theodoridis and Koutroumbas 2008). A brief overview of this technique is given below for the benefit of the reader. Given an l × n matrix X of rank r (r ≤ min{l, n}), there exist unitary matrices U and V of dimensions l × l and n × n respectively, such that

$$\begin{equation} X = U\left[ {\begin{array}{@{}*{20}c@{}} {\Lambda ^{1/2} } & {\rm O} \\ {\rm O} & 0 \\ \end{array}} \right]V^T , \end{equation} \tag{ 10 }$$

where Λ^1/2 is the diagonal matrix with λ_i^1/2 elements, and λ_i are the r nonzero eigenvalues of the associated matrix X^HX. By means of these matrices, it is possible to expand an approximate representation of X using less than r terms in reconstruction. The error in the approximation is the sum of the unused singular values (Theodoridis and Koutroumbas 2008). Unlike principal component analysis (PCA), the SVD components may show some degree of correlation. To decompose the features matrix (X), we used all records from healthy and diabetic subjects together. The size of the feature matrix for each group being 8 × 23, by horizontally concatenating these matrices a size of 8 × 46 is obtained for X. The sum of the first three singular values was found to be more than 90.1% (for healthy subjects) and 99.2% (for diabetics) of the total sum for the eight singular values. In order to determine the optimum number of SVD components, this number was varied and it was observed that using four SVD components instead of three did not improve the performance of the classification. Further increasing the number of components to above four resulted in even worse performance. It was concluded that the optimum number of SVD components in this feature space would be three, and the eight primary features from (9) were reduced to three. To reduce the number of features from eight to three, the classical SVD method was employed using the transformation matrix T as follows:

$$\begin{equation*} T = \left[ {\begin{array}{@{}ccccccccc@{}} { - 7} & {559} & {56} & { - 678} & { - 99} & {401} & {183} & { - 144} \\ { - 11} & { - 715} & {29} & { - 116} & { - 233} & {511} & {247} & { - 314} \\ { - 4} & {412} & { - 31} & {652} & { - 286} & {214} & {149} & { - 504} \\ \end{array}} \right] \times 10^{ - 3}. \end{equation*}$$

The transformation matrix T was computed once and then used to generate feature vectors with lower dimension. The product of transformation matrix T and the 8 × 46 features matrix X produced the new features matrix with a size of 3 × 46, reducing the number of features to three for each subject. Interestingly enough, a 3D plot of these SVD-generated features (figure 5) revealed that each group's features are contained within a plane.

Zoom In Zoom Out Reset image size

Using the least-squares method, the equation for the plane representing the locus of the SVD-reduced parameters in healthy subjects was obtained as

$$\begin{equation} z = - 1.7 - 0.28x - 0.84y\end{equation} \tag{ 11 }$$

with a sum of squared errors (SSE) of 0.21, whereas for the diabetic group the fitted plane is

$$\begin{equation} z = 1.48 + 0.25x - 0.58y\end{equation} \tag{ 12 }$$

with an SSE of 0.04.

2.8. Statistical analysis

2.9. Bayesian classifier for two classes with normal distribution

Figure 6 shows the mean of the first three SVD-generated features and their 95% confidence intervals. A p-value <0.05 was obtained for the third SVD component, whereas the second component resulted in p <0.001. Table 1 summarizes the numerical results.

Zoom In Zoom Out Reset image size

Next the accuracy, sensitivity and specificity were computed according to the following equations:

$$\begin{equation} {\rm Accuracy}{\rm (performance)} = \frac{{{\rm TP} + {\rm TN}}}{{{\rm TP} + {\rm FP} + {\rm FN} + {\rm TN}}} \end{equation} \tag{ 13 }$$

$$\begin{equation} {\rm Sensitivity} = \frac{{\rm TP}}{{{\rm TP} + {\rm FN}}} \end{equation} \tag{ 14 }$$

$$\begin{equation} {\rm Specificity} = \frac{{\rm TN}}{{{\rm TN} + {\rm FP}}}, \end{equation} \tag{ 15 }$$

where true positive (TP) and false positive (FP) are the number of correctly and incorrectly classified patients, respectively. True negative (TN) and false negative (FN) are the number of correctly and incorrectly classified healthy subjects, respectively.

Table 2 shows the results when a naive Bayes classifier was applied on the first three significant components of SVD (SVD1, SVD2 and SVD3) for all subjects. When subjects are separated according to their gender (second and third rows in table 2), a considerable improvement is realized in the classification. It can be noted that the overall accuracy in the male group is higher than that for the female group.

The main objective of this research was to investigate changes in the upper limb vascular system following a stimulus to endothelial layer cells by induced shear stress. The shear stress was created by provoking local ischemia, and the effect of the resulting FMD was investigated amongst 23 diabetics (with at least one cardiovascular risk factor) and 23 healthy subjects. This evaluation was performed using a differential modeling of cardiovascular system that uses bilateral PPG. Results point to the fact that when the vascular system is stimulated, by measurement of its response, its state can be efficiently identified with a relatively good accuracy (93.5%). When gender is taken into consideration, results are further improved.

Although ultrasound FMD is a relatively established method to assess EF in research settings, there are still many obstacles rendering its clinical use unacceptable for physicians. A technique that could potentially eliminate the need for a very costly and operator-dependent ultrasound system can be considered a step forward to establish FMD as a clinical test procedure. In spite of many research studies toward this replacement (Zahedi et al 2008, Selvaraj et al 2009), none has yet achieved high correlation (e.g. above 80%) between PPG-FMD and US-FMD. Generally, in these studies, the PPG waveform is morphologically investigated to estimate the BA diameter. One hypothesis regarding the major hindrance to achieve this objective is the auto-regulation process in small arteries at the fingertip (Klabunde 2012) that causes the amplitude of the PPG response to be tens orders of magnitude greater than BA diameter changes (Zahedi et al 2008). In the present work, the PPG-FMD process was modeled from a systemic point, and subsystems contributing to shape blood volume changes at the fingertip were part of the whole model. Results show that the effect of the auto-regulation on the PPG during FMD seems interestingly taken into account by our model, and corresponding classification results confirm the hypothesis that the auto-regulation mechanism is heavily affected by diabetes.

Our work shows some improvement over a previous tentative to investigate the effect of diabetes on PPG (Shariati et al 2008) using a differential vasculature model. This highlights the importance of using a stimulus to elicit a response allowing for a discrimination between the two groups of subjects. Another interesting point is that whereas gender-based classification resulted in better classification in our work, in Shariati et al (2008) gender-based classification results in lower accuracy.

One of the limitations of this work that needs to be overcome in future investigations is to use age-matched groups. Indeed, age has been shown to play a relatively important role in the genesis of vascular dysfunction (Lakatta 1990) and needs to be eliminated as a confounding factor by selecting age-matched subjects. In practice, this means having access to a very large pool of volunteers and choosing those who fit the age criteria.

There are many ways to improve feature discrimination between classes, namely feature selection and generation. While PCA is a feature generation method based on second-order statistics, SVD employs a limited number of observations to construct a transformation matrix. This may explain why using PCA leads to poor results. By means of SVD, over 90% of the total sum of singular values is compacted into just three features and this may explain our high classification results.

Most of the subjects in this study are at early stages of the disease (duration of diabetes is less than 2 years), and our results confirm the hypothesis that the changes in cardiovascular system could be detected by our differential modeling. Although here we concentrate only on diabetic subjects, other cardiovascular risk factors could be considered as well to investigate the detectability of changes on the EF at early stages.

The dysfunction of endothelial cells is a sign of cardiovascular disease that precedes other clinical symptoms. Thus, the investigation of this response helps to diagnose these diseases at early stages. Our results confirm the hypothesis that cardiovascular impairments due to diabetes may be detected through a non-invasive and inexpensive PPG-FMD test, using differential modeling of the upper limb cardiovascular system.

The high sensitivity (100%) of our model when gender is taken into account warrants more investigation by especially increasing the sample size. However, this method looks promising in the detection of endothelial dysfunction.

Further work is needed to investigate the accuracy of this method in larger sample sizes. Besides diabetic subjects, a future step would be to investigate this model on subjects with other cardiovascular risk factors. Finally, due to the effect of age, age matching should be ensured to further ascertain the validity of the results.

The authors would like to thank the reviewers for sharing their comments, which considerably improved the quality of the paper. They are also grateful to the National University of Malaysia (Universiti Kebangsaan Malaysia), Faculty of Engineering and Built Environment for kindly providing access to the data for this study. This work has been partially funded by the Iranian National Science Foundation (INSF) grant number 87041735.