Assessment of cardiovascular function from multi-Gaussian fitting of a finger photoplethysmogram

The susceptibility of the PPG signal to ambient light at the photo-detector, poor blood perfusion of the peripheral tissues and motion artifacts (Sukor et al 2011) is well known. In order to minimize the influence of these factors in the subsequent phases of the PPG analysis, a pre-processing stage is required. At this stage, high frequency components in the acquired signal, such as the interference of the power sources, are removed using a low-pass Butterworth filter with a 18 Hz cut-off frequency (Chan et al 2007), ensuring the preservation of the physiological relevant information (usually bellow 15 Hz (Kamal et al 1989)). Additionally, low frequency components, commonly associated with changes in capillary density and venous blood volume (Reisner et al 2008), temperature variations, acquisition instruments or subject movements (Pilt et al 2008), were suppressed by subtracting a two-second window moving average filtered version to the original PPG signal (Chan et al 2007).

The main objective in this step is to segment a PPG signal stream into individual PPG pulses per heart beat. From each detected PPG pulse, the morphological characteristics are interred afterwards.

To detect the PPG pulses, an approach similar to Sun et al (2005) was applied. Firstly, the PPG signal is differentiated using a five-point digital differentiator (Abramowitz and Stegun 2012) (applying equations (1)–(4) in parallel to the PPG signal), resulting in first- to fourth-order derivatives (d1_ppg, d2_ppg, d3_ppg and d4_ppg).

$$\begin{eqnarray}&&d{{1}_{\text{ppg}}}={{f}^{'}}(t)=\frac{f(t-2h)-8f(t-h)+8f(t+h)-f(t+2h)}{12{{h}^{2}}}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&d{{2}_{\text{ppg}}}={{f}^{\prime \prime}}(t)=\frac{-f(t-2h)+16f(t-h)-30f(t)+16f(t+h)-f(t+2h)}{12{{h}^{2}}}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&d{{3}_{\text{ppg}}}={{f}^{\prime \prime \prime}}(t)=\frac{-f(t-2h)+2f(t-h)-2f(t+h)+f(t+2h)}{2{{h}^{3}}}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&d{{4}_{\text{ppg}}}={{f}^{\prime \prime \prime \prime}}(t)=\frac{f(t-2h)-4f(t-h)+6f(t)-4f(t+h)+f(t+2h)}{{{h}^{4}}}\end{eqnarray} \tag{ 4 }$$

where f is the PPG time series, t is the time index and h is the sampling time.

From the analysis of the d1_ppg, the local maxima (d1_ppg_lmax) with absolute amplitude greater than a threshold ThR are detected. Here, the ThR is selected based on an adaptive thresholding of the d1_ppg data cumulative histogram (calculated within 10 s time windows), and was defined as the greater value bellow which 90% of the observations are found.

Following the work proposed by Chan et al (2007), the PPG pulse onset was considered to be the point d1_ppg presenting a rapid inflection (just before d1_ppg_lmax), which corresponds to a maximum in the d3_ppg (d3_ppg_lmax). The motivation behind this approach was firstly reported by Cook (Cook 2001, Wisely and Cook 2001), who observed a significant resemblance between the PPG signal's first derivative and the arterial flow waveform.

The detection of the d3_ppg_lmax characteristic point was accomplished in two phases: (1) detection of the d3_ppg local minima (d3_ppg_lmin) corresponding to the d1_ppg local maxima (see figure 2) and (2) identification of the peak with greater amplitude (d3_ppg_lmax) prior to the previously identified most relevant valley.

Zoom In Zoom Out Reset image size

The first step of the modeling process is an additional preprocessing of the PPG pulses, which may still exhibit minor baseline and amplitude fluctuations and can negatively affect the subsequent modeling phases. In order to reduce these potential sources of error, the linear trend is removed and the amplitude of each PPG pulse is normalized to the unit. Thus, we aim to reduce the variability of the model parameters resulting from exogenous factors, improve the optimization procedure and avoid the normalization of the model parameters in a latter stage.

In the second step, the PPG pulse is separated into systolic and diastolic phases. The systolic phase, primarily associated with the ventricular ejection, was defined as the portion between the onset of the PPG pulse and the dicrotic notch. The diastolic phase, resulting from the reflections and re-reflections of the main PPG pulse in the systemic vascular network, was defined as the portion of the PPG pulse between the dicrotic notch and the pulse offset. The two phases of the PPG pulse are illustrated in figure 3.

Zoom In Zoom Out Reset image size

The detection of the reference point correspondent to the dicrotic notch was defined as local maxima of the d2_ppg (d2_ppg_lmax) within a [0.2–0.4]s interval and its boundaries as the negative-to-positive/positive-to-negative zero crossings of the d2_ppg (Rubins 2008). Although this strategy can be effective in PPG pulse morphologies where the dicrotic notch is highly prominent, in cases where the pulse is severely damped this task can be deeply compromised due to the existence of multiple d2_ppg_lmax and the absence of zero crossings in an acceptable temporal window (i.e. very far from d2_ppg_lmax—an example in shown in figure 4). In such cases, the identification of the dicrotic notch boundaries relies on the identification of the points in d2_ppg with the steepest slope, i.e. the negative-to-positive/positive-to-negative in the d4_ppg closer to the reference point. The selection of the dicrotic notch among multiple candidates is performed based on a joint analysis of the candidate characteristics, i.e. their amplitude, area and length.

Zoom In Zoom Out Reset image size

In our previous work (Couceiro et al 2012), we modeled the PPG pulse as a mixture of four Gaussian functions, related with the main pulse and corresponding reflections in the arterial path. However, depending on the morphology of the PPG pulse, this model may lead to inaccurate approximations, especially in rapid inflections close to the PPG pulse onset during the systolic rise. Based on these observations and aiming at reduction of the error between the proposed PPG model and the data, the main forward pulse is modeled as a mixture of two Gaussian functions instead of just one. Thus, by increasing the number of degrees of freedom of the proposed model, we expect to decrease the error between the model and the data and to determine the PPG pulse components more accurately. The formulation of the proposed five-Gaussian model is presented as follows:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{f}_{m}}\left(t,{{\beta}_{j}}\right)=\underset{j=1}{\overset{5}{\sum}}\,{{g}_{j}}\left(t,{{\beta}_{j}}\right),\text{}\text{ for}\text{}{{\beta}_{j}}={{\left\{a,b,c\right\}}_{j}} \\ {{g}_{j}}\left(t,{{\beta}_{j}}\right)={{a}_{j}}\text{}\text{}{{\text{e}}^{\frac{-{{\left(t-{{b}_{j}}\right)}^{2}}}{2{{c}_{j}}^{2}}}} \\ \end{array}\end{eqnarray} \tag{ 5 }$$

where the parameters a_j, b_j and c_j correspond to the amplitude, location and width of the Gaussian function j. Here, the first two Gaussian functions (${{g}_{1\_2}}={{g}_{1}}+{{g}_{2}}$) correspond to the ventricular ejection and the third Gaussian function (g₃) is related to the first pulse reflection at the junction between the thoracic and abdominal aorta, presented in figure 5. The fourth (${{g}_{4}}$) and fifth (${{g}_{5}}$) Gaussian functions derive from the reflection at the juncture between abdominal aorta and common iliac arteries (Baruch et al 2011) and minor reflections and re-reflections in the systemic structure, respectively (see figure 5).

Zoom In Zoom Out Reset image size

According to Wang et al (2013), the positions of the a, b, c, d and e waves (here called W_a, W_b, W_c, W_d and W_e waves and illustrated in figure 6) defined by DDA (Takazawa et al 1998) can provide an insight about the locations of the individual components of the PPG pulse. Additionally, from the inspection of the PPG pulse morphology and the corresponding second derivative, we observed that the negative peak after W_e (i.e. W_f) can provide information about the location of the second reflection wave and, therefore, it was included in the present work. Based on the waves identified figure 6, the initial parameters and corresponding boundaries of the MG model were determined using a Heuristic approach similar to the one proposed by Wang et al (2013). The initial parameters and boundaries of the proposed model are presented in table 1.

Zoom In Zoom Out Reset image size

To adjust the model parameters a least squares method combined with the interior point algorithm (Waltz et al 2006) has been adopted. Let $f(n),~n\in \left[1,N\right]~$ be the normalized PPG pulse, where N is the total number of samples composing $f(n)$. Using the least squares method, the goal is to find the solution, that minimize the sum of the squared residuals as presented bellow:

$$\begin{eqnarray}&&{{f}_{\text{MSE}}}(\beta )=\frac{1}{N}\underset{n=1}{\overset{N}{\sum}}\,{{\left[f(n)-{{f}_{m}}(n,\beta )\right]}^{2}}\end{eqnarray} \tag{ 6 }$$

where ${{f}_{m}}(n,\beta )$ is the MG model with the set of parameters ${{\beta}_{j}}={{\left\{a,b,c\right\}}_{j}},~j=1,\ldots ,5$. Using the interior point algorithm (Waltz et al 2006), we aim to solve following nonlinear constrained optimization problem

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \text{min}\,{{f}_{\text{MSE}}}(\beta ) \\ {} \\ \text{subject}\text{ to}\,\,g(\beta )\leqslant 0\text{and}\text{lb}\leqslant \beta\leqslant \text{ub}\,\, \\ {} \\ \end{array}\end{eqnarray} \tag{ 7 }$$

where $g(\beta )\leqslant 0$ and $lb\leqslant \beta \leqslant ub$ are the inequality constrains and the boundaries to which the parameters are subject. By defining the parameters initial values we ensure that the optimization starting point is closest to the optimal solution (based on a priori information). Additionally, constraining the parameters to pre-specified boundaries and relationships between them, we aim to limit the number possible solutions and assure the physiological foundations and interpretability of the model, as well as to assure its robustness and stability. Based on the study proposed by Baruch et al (2011), the inequality constraints are summarized bellow:

• ${{a}_{1}}-{{a}_{2}}\leqslant 0\,\wedge \,{{a}_{1}}-{{a}_{3}}\leqslant 0$: the main objective of ${{g}_{1}}$ is to increase the fitness of the proposed model in the first portion of the systolic rise. Since it belongs to the main forward wave, we confine the amplitude of ${{g}_{1}}$ must be smaller than the following two (${{g}_{2}}$ and ${{g}_{3}}$).
• ${{a}_{3}}-{{a}_{2}}\leqslant 0\text{}\text{}\wedge {{a}_{4}}-{{a}_{2}}\leqslant 0\wedge {{a}_{5}}-{{a}_{2}}\leqslant 0$: the forward pulse has the highest amplitude (Baruch et al 2011) and therefore the amplitude of ${{g}_{2}}$ must be greater than ${{g}_{3}}$, ${{g}_{4}}$ and ${{g}_{5}}$.
• ${{a}_{5}}-{{a}_{4}}\leqslant 0$: since ${{g}_{5}}$ is a result of reflections and re-reflections in the arterial pathway, this wave is expected to have very low amplitude, i.e. the smallest within the diastolic phase. Therefore we restrict the amplitude of the last reflected wave (${{g}_{5}}$) smaller than the second reflected pulse (${{g}_{4}}$).
• ${{b}_{1}}\leqslant {{b}_{2}}<{{b}_{3}}\leqslant {{b}_{4}}\leqslant {{b}_{5}}$: in order to guarantee that the aforementioned relationships prevail, the location of the Gaussian functions must be increasingly larger.

As the main pressure wave, created by the ejection of blood from the left ventricle, travels in the arterial pathway, it suffers from dampening, which affects both the rise and decay of the pulse. Additionally, the reflection of the main pulse in the two major reflection sites also causes drastic changes on the diastolic phase of the measured PPG pulse (Chan et al 2007). In order to avoid these potential sources of error, one implemented a data driven approach that resorts on the analysis of the components corresponding to the systolic phase of the PPG pulse, herein called as the systolic model (${{g}_{\text{sys}}}$). ${{g}_{\text{sys}}}$ was defined as the sum of the Gaussian functions corresponding to the main forward wave (${{g}_{1}}$ and ${{g}_{2}}$) and late systolic peak (${{g}_{3}}$). Since the morphological characteristics of the pulse wave (at the beginning of the aorta) change as the pulse travels the pathway, it is expectable that the time span between the onset and offset of the ${{g}_{\text{sys}}}$ model do not exactly match the time span between the onset and offset of the left ventricular ejection. Therefore, we aim to investigate which are the characteristic points at the beginning and end of the ${{g}_{\text{sys}}}$ model that better represent the LVET. For that purpose, a data driven approach based on a multi-derivative analysis of the ${{g}_{\text{sys}}}$ model has been implemented. ${{g}_{\text{sys}}}$ was differentiated using (1)–(3), resulting in the derivatives up to the 3rd order and six LVET estimates (LVETc_i) from the analysis of the systolic model derivatives:

• LVETc₁—time span between systolic (D3_sp) and diastolic (D3_dp) peaks of the systolic model third derivative.
• LVETc₂—time span between systolic peak of the systolic model third derivative (D3_sp) and the notch in the systolic model 1st derivative (D1_nt).
• LVETc₃—time span between systolic (D2_sp) and diastolic (D2_dp) peaks of the systolic model second derivative.
• LVETc₄—time span between systolic peak of the systolic model second derivative (D2_sp) and diastolic peak of the systolic model third derivative (D3_dp).
• LVETc₅—time span between systolic peak of the systolic model third derivative (D3_sp) and diastolic peak of the systolic model second derivative (D2_dp).
• LVETc₆—time span between systolic peak of the systolic model second derivative (D2_sp) and the notch in the systolic model first derivative (D1_nt).

The characteristic points and respective LVET estimates are illustrated in figure 7.

Zoom In Zoom Out Reset image size

Several parameters have been suggested in the literature as potential surrogates of blood pressure and vascular tone changes (O'Rourke et al 2001, Millasseau et al 2002, Allen 2007, Baruch et al 2011, Meyer et al 2011, Muehlsteff et al 2012, 2013). In the current paper, we investigate two major parameters called stiffness index (SI) and reflection index (RI). The SI is associated with the velocity of a pulse wave in large arteries (DeLoach and Townsend 2008) and large artery stiffness (Millasseau et al 2002). It also correlates with pulse pressure (Baruch et al 2011) and is defined by the time span between the forward (P_s) and reflected waves (P_d), as illustrated in figure 8, that is:

$$\begin{eqnarray}&&\text{SI}={{B}_{d}}-{{B}_{s}}\end{eqnarray} \tag{ 8 }$$

Zoom In Zoom Out Reset image size

Additionally, the RI, associated with small artery stiffness (DeLoach and Townsend 2008) and changes in vascular tone (Chowienczyk et al 1999, Laucevicius et al 2002), is defined by the ratio between the amplitudes of both waves the forward (P_s) and reflected waves (P_d), described as follows:

$$\begin{eqnarray}&&\text{RI}=\frac{{{A}_{d}}}{{{A}_{s}}}\end{eqnarray} \tag{ 9 }$$

Additionally, we also aim to investigate the temporal and amplitude relationship between the main forward wave (P₁) with the late systolic wave (P₂) and the reflected waves (P_d) described by the four additional indexes (T_{1_2}, T_{1_d}, R_{1_2}, R_{1_d}), presented from (10)–(13).

$$\begin{eqnarray}&&{{T}_{1\_2}}={{B}_{2}}-{{B}_{1}}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{{T}_{1\_d}}={{B}_{d}}-{{B}_{1}}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{{R}_{1\_2}}=\frac{{{A}_{2}}}{{{A}_{1}}}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{{R}_{1\_d}}=\frac{{{A}_{d}}}{{{A}_{1}}}\end{eqnarray} \tag{ 13 }$$

While the indexes T_{1_d}, and R_{1_d}, are closely related with the SI and RI, the indexes R_{1_2} have been associated to changes the impedance mismatch of the aortic artery and is closely related with augmentation index, which has been associated aging and arterial stiffening (O'Rourke et al 2001, Allen 2007, Baruch et al 2011).

Two databases were used in the current work. The evaluation of the LVET detection algorithm was performed on a database collected at the Centro Hospitalar de Coimbra (CHC), while for the investigation of the SI, T_{1_d}, T_{1_2}, RI, R_{1_d} and R_{1_2}, a database provided by the Heinrich-Heine University Hospital, Düsseldorf was used. The protocols used for the collection of these databases are outlined bellow.

3.1.1. Protocol for the assessment of cardiac function surrogates.

The data considered in the present paper for the assessment of cardiac function was collected in a study conducted in CHC aiming at the simultaneous acquisition of PPG and ECHO.

In this study, patients of two distinct groups were enrolled: (1) 33 healthy subjects and (2) 35 subjects suffering from various cardiovascular diseases (CVDs) including coronary artery disease, heart failure, and comorbidities such as hypertension. All volunteers gave written informed consent to participate in this study, which was authorized by the CHC's ethical committee in 2010 under the protocol 'Assessment of cardiac function using heart sounds, ICG and PPG'. The characteristics of the population considered in the present study are summarized in table 2.

Table 2. Characteristics of the volunteers enrolled in the study for the assessment of cardiac function.

	Healthy (#33)	CVD (#35)
Age	29   ±   8	58   ±   17
Weight	74.3   ±   11.1	70.6   ±   12.6
BMI	24   ±   2	25   ±   3
Male/female	32/1	20/15

The ECG and PPG signals were recorded with an HP-CMS monitor and stored on a laptop at a sampling frequency of 500 Hz and 125 Hz, respectively. ECHO (in Doppler mode, with ECG) was conducted using a Vivid 7 system from General Electric, which produces images with a temporal resolution of 500 Hz (figure 9). The ECG signals acquired simultaneously by the two systems described above were used to synchronize the signals streams.

Zoom In Zoom Out Reset image size

The measurement protocol was supervised and conducted by an authorized medical specialist. It consisted of several acquisitions of ECHO in Doppler mode and PPG collected at the right-hand index finger using an infrared transmission finger probe.

A clinical expert annotated the ECHO-doppler images. The time instant corresponding to the opening of the aortic valve was defined as onset of the left ventricular outflow lobe, while the closure of the aortic valve was defined immediately before the closing click onset, produced by the residual reflux after the closure of the aortic valve cusps. An example is shown in figure 9.

In this study a total of 2081 heartbeats have been annotated, of which 1109 beats correspond to healthy volunteers and 972 beats correspond to CVD subjects.

3.1.2. Protocol for the assessment of blood pressure and vascular tone surrogates.

The synchronization of the data obtained during the before mentioned studies was performed using an algorithm based on the minimization of the error between the RR intervals time series assessed from the ECGs of both systems, using the minimum least squares error criterion. Let RR₁(k), k = 1,..., K, and RR RR₂(m), m = 1,..., M (K  <  M) be the RR intervals time series of the ECGs to be synchronized. The time instant t at which both time series are synchronized is:

$$\begin{eqnarray}&&t=\underset{w=0,\ldots ,M-K}{{\text{argmin}}}\,{{f}_{\text{sync}}}(w)\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{{f}_{\text{sync}}}(w)=\frac{1}{K}\underset{k=1}{\overset{K}{\sum}}\,~{{\left(RR1~(k+w)~RR2~(k)\right)}^{2}}\end{eqnarray} \tag{ 15 }$$

where ${{f}_{\text{sync}}}$ is the synchronization time series.

The algorithms' performance for LVET inference was analyzed for the whole dataset described in section 3.1.1, referred to as global dataset in the following, and in the two subsets that compose it. These are the subset of healthy volunteers, called the healthy subset, and the subset of volunteers with cardiovascular diseases, called the CVD subset.

To evaluate the performance achieved by the proposed algorithm, four performance metrics have been adopted. Let LVET_MEAS:{LVETc₁,..., LVETc₆, LVET_CHAN} be the measured beat-to-beat LVET values, corresponding to estimate 1–6 and the LVET obtained using the algorithm proposed by Chan et al (2007), respectively, and LVET_ECHO be the reference beat-to-beat LVET values assessed using ECHO-doppler. The abbreviation 'Error' stands for the error between measured and reference values ($\text{LVE}{{\text{T}}_{\text{MEAS}}}-\text{LVE}{{\text{T}}_{\text{ECHO}}}$), while 'Abs. Error' concerns to the absolute estimation error ($|\text{LVE}{{\text{T}}_{\text{MEAS}}}-\text{LVE}{{\text{T}}_{\text{ECHO}}}|$). The 'Abs. Error Perc.' stands for the percentage of absolute estimation error ($|\text{LVE}{{\text{T}}_{\text{MEAS}}}-\text{LVE}{{\text{T}}_{\text{ECHO}}}|/\,\overline{\text{LVE}{{\text{T}}_{\text{ECHO}}}}$). The agreement between measured LVET and LVET_ECHO was analyzed by assessing the correlation coefficients between the two variables in each dataset (Global, healthy and CVD). Additionally, the correlation coefficients were also calculated for each volunteer and the average and standard deviation over each dataset was calculated. The Pearson correlation coefficient was calculated if the variables were normally distributed, while Spearman's correlation coefficients were used otherwise. The p-values corresponding to the coefficients shown in table 4 were low (p  <  0.05), enabling the rejection of the hypothesis of no correlation.

The agreement between $~\text{LVE}{{\text{T}}_{\text{MEAS}}}$ and $\text{LVE}{{\text{T}}_{\text{ECHO}}}$ was also analyzed using regression plots (figure 10) and a Bland–Altman plot (figure 11). Error distributions were tested for Gaussianity by using the Kolmogorov–Smirnov test. Accordingly, statistical analysis was performed using the paired Student test and the two-sided Wilcoxon signed rank test and the values shown in table 4 have been found to be statistically relevant (p  <  0.001).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In table 4 we present the results achieved by our algorithm approach for all extracted LVET estimates. One observes that the best estimation performance was achieved by the LVETc₆ estimate (measured between D2_sp and D1_nt), with an absolute estimation error of 15.41   ±   13.66 ms. This error translates to a percentage of error of 5.53   ±   4.9% which is comparable to the inter-observer variability (3.6%–5.4%) of Doppler echocardiographic LVET measures reported in literature (Sohn and Kim 2001, Tekten et al 2003, Colan et al 2012). The performance of the remaining estimates in the global dataset, are much less good ranging from 26.13   ±   14.76 to 75.96   ±   23.78 ms. An exception is observed for the LVETc₄ estimate (measure between D2_sp and D3_dp) that obtained an absolute estimation error of 18.72   ±   16.51 ms, representing an increase in the estimation error percentage of about 1% (3 ms). The correlation between the measured LVET estimates and the echocardiographic LVET reference was high and ranged from 0.68 (LVETc₅) to 0.78 (LVETc₆).

The results regarding the evaluation of the proposed LVET estimates in both healthy and CVD subsets are also presented in table 4. In the healthy subset, we observed a decrease in the absolute estimation performance of about 5 ms for both LVETc₄ and LVETc₆, representing a reduction of about 1–2%. Contrarily, a decrease in the performance of the remaining estimates was observed for this subset, with an increase in the estimation error ranging from 3 to 5%. The agreement between both measure and reference LVET's increased for all the LVET estimates, with an exception to LVETc₆, where a decrease of about 6% was observed. Contrarily to the results observed in the healthy subset, in the CVD subset we observed a decrease in the performance achieved by LVETc₄ (23.99   ±   17.67 ms) and LVETc₆ (20.63   ±   15.75 ms) and a decrease in the absolute estimation error of remaining LVET estimates. Moreover, an increase in the agreement between measured and reference LVET values was also observed in this study group, being LVETc₆ the estimate with highest correlation coefficient (ρ = 0.85).

The correlation coefficients calculated in a patient-by-patient basis presented results similar (but lower,≈0.15 less) to those calculated in each dataset, showing a lower agreement between the measured LVET values and the LVET_ECHO, which is mainly explained by the absence of marked trend within each volunteer measured values.

Comparing the best LVET estimate (LVETc₆) with the LVET obtained using the algorithm proposed by Chan et al (2007), one observes that the best estimation performance was achieved by our method, which obtained better estimation performance (15.41   ±   13.66 ms versus 23.01   ±   14.60 ms) and higher correlation (0.78 versus 0.75). In the healthy context, a significant improvement as been achieved by LVETc₆ (Abs. Error: 10.88   ±   8.79 ms and ρ = 0.72) compared to the algorithm by Chan et al (2007) (Abs. Error: 31.03   ±   12.66 m and ρ = 0.71). In contrast, Chan's method performs better for CVD patients compared to our implementation (13.48   ±   10.68 ms versus 20.63   ±   15.75 ms). Finally, Chan's algorithm achieved a better correlation coefficient for CVD patients (ρ = 0.88 versus ρ = 0.85). Although our algorithm achieves higher estimation error and lower correlation in the CVD dataset, the results achieved in a balanced population (global dataset) with the similar number of beats from healthy and CVD subjects suggests that our algorithm is able to perform better across large populations.

In figure 10 we present the regression plots corresponding to the best LVET estimate (LVETc₆) in both healthy and CVD subsets. For the healthy subset (figure 10(a)), the best linear fit corresponds to a 0.94 slope, which is very close to unity. This agrees with the low estimation error achieved by LVETc₆ and the high correlation value presented in table 4. For the CVD subset figure 10(b)), a similar slope of 0.94 was achieved, but with a greater dispersion of the reference/measure LVET pairs.

The agreement between the LVETc₆ and LVET_ECHO was also assessed via Bland–Altman plot, presented in figure 11 for the global dataset. The horizontal solid line represents the mean error while the dashed lines denote the level of agreement between LVETc₆ and LVET_ECHO, i.e. mean  ±  standard deviation and mean  ±2 standard deviation. One can observe that the average error is close to zero and the estimation error is not evenly distributed around the whole range of LVET values, particularly in the higher range where a greater dispersion is shown.

The presented results show that the characteristic points that better define the onset and end of the ventricular ejection on the PPG pulse are the systolic peak defined in the second derivative and the notch defined the first derivative. Unsurprisingly, these two points do not match exactly the onset and offset of the PPG main forward wave resultant from the ventricular ejection. This can be explained by the temporal and morphological changes experienced by the main pulse as it travels from the left ventricle to the peripheral sites. The reduction of arterial compliance (O'Rourke 1995) and the increased peripheral resistance (Allen 2007), due to the changes in the arterial wall composition and the tapering of the arteries, dampens the main pulse, and therefore a temporal shift is observed in the onset and offset of the PPG pulse. Instead, the characteristic points corresponding to the boundaries of the ejection in the PPG systolic model are determined by the characteristic points corresponding to the maximum acceleration at the systolic rise of the systolic model and by the minimum velocity during its decay. Contrarily, the remaining LVET estimates calculated form the characteristic points closer to the onset/offset of the systolic model present a higher error, resulting from an overestimation.

Our results also indicate a discrepancy in the estimation of LVET within healthy and CVD patients. In the healthy dataset, the absolute estimation error achieved by LVETc₆ was very low (ap. 10.9 ms), while in the CVD the absolute error almost doubled (ap. 20.6 ms) accentuated by the underestimation of the LVET (error:  −17.65   ±   19.04 ms). One possible reason for the underestimation of LVET in the CVD population is the dependence of the pulse transmission on the PWV. One common characteristic observed in aging and in populations with cardiovascular diseases, is the increase of the pulse wave velocity, mainly due to the increase of the vascular resistance and stiffness, which leads to a faster transmission of reflected waves along the arterial pathway and early appearance in the PPG pulse (O'Rourke et al 2001, Allen 2007). With the aging of the arteries (and the presence of atherosclerosis), there is a progressively return of the reflected waves, which can arrive so early during systole that it becomes difficult to distinguish the two phases (Millasseau et al 2002). This phenomenon is believed to negatively affect the proposed algorithm in two aspects. The first is the increase in the uncertainty in the identification of the systolic and diastolic phase boundaries, which consequently affects the remaining modeling process. The second is the overestimation of the Gaussian parameters related with the reflected waves in the diastolic phase, due to the attempt of the model to fit the increased (in amplitude and length) falling edge of the diastolic phase. As a consequence of this overestimation, the amplitude and length systolic model is diminished leading to an underestimation of the LVETc₆. Another characteristic that is believed to accentuate the problematic is the elongation of the systolic rising edge, commonly associated with aging. As a result of this elongation, the characteristic points corresponding to the maximum acceleration of the systolic model appear at a later instant and therefore a decrease in the LVET estimate is observed.

In this section, we analyze the relationship between the extracted blood pressure and vascular tone surrogates with the reference blood pressure and total peripheral resistance values. The relationship was evaluated in the dataset composed by all volunteers (global dataset) and in the subset composed by volunteers where greater blood pressure changes are observed, i.e. the tilt positive group (tilt positive dataset). The relationship between variables was evaluated for the upright position phase of the HUTT protocol. The extracted reference values used in the present analysis are: systolic (SBP), diastolic (DBP) and mean blood pressure (MBP), pulse pressure (PP) and total peripheral resistance index (TPRI). These parameters were compared with the extracted surrogates, i.e. SI, T_{1_d}, T_{1_2}, RI, R_{1_d} and R_{1_2}, in order to identify the parameters that better characterize blood pressure and vascular tone changes, and their level of agreement with those parameters.

Let m = 1,..., M with M = 43 be the volunteers, x_meas = {SI, T_{1_2}, T_{1_d}, RI, R_{1_2}, R_{1_d}} the measured parameters and x_ref = {SBP, DBP, MBP, PP, TPRI} the reference parameters, assessed in a beat-to-beat basis. The agreement (${{\rho}_{p}}$) between the measured and reference parameters for each volunteer using the Pearson correlation coefficient was defined as:

$$\begin{eqnarray}&&{{\rho}_{\text{p}}}\left({{x}_{\text{meas}}},{{x}_{\text{ref}}}\right)=\frac{1}{M}\underset{m=1}{\overset{M}{\sum}}\,|\rho \left({{x}_{\text{meas}}},{{x}_{\text{ref}}}\right)|\end{eqnarray} \tag{ 16 }$$

In tables 5 and 6 we present the Pearson correlation coefficients ρ_p achieved by each extracted features and the reference parameters, for the global and tilt positive datasets. The p-values corresponding to the coefficients shown in tables 5 and 6 were low (p  <  0.05), enabling the rejection of the hypothesis of no correlation.

From tables 5 and 6, one can observe that the parameter achieving the best correlation coefficient in the global dataset is the T_{1_2}, which presented a correlation with TPRI, SBP of 0.45 and 0.43, respectively. Additionally, the parameter SI and T_{1_d} were also moderately correlated with the reference blood pressure measurements and with the TPRI, with values ranging from 0.37 (between T_{1_d} and TPRI) to 0.42 (between SI and DBP). Similarly, the parameters RI and R_{1_d} also showed a good agreement with the reference parameters (excluding PP) being the highest correlation with TPRI (ρ_p = 0.43 and ρ_p = 0.42, respectively). The parameter R_{1_2} showed to be weakly associated with each of the reference parameters, with ρ_p ranging from 0.21 to 0.26.

Regarding the results in the tilt positive group, one can observe a marked increase in the agreement between the majority of the measured and reference parameters. In this group of patients, the highest agreement was found between the SI and MBP (ρ_p = 0.57), SBP (ρ_p = 0.57) and DBP (ρ_p = 0.54). Not surprisingly, the parameter T_{1_d} also presented a high agreement with MBP and SBP (ρ_p = 0.53). The parameter most closely associated with TPRI in the tilt positive volunteers was the SI (ρ_p = 0.51). The parameter T_{1_2} showed to be more closely associated with the SBP, PP and TPRI, with a level of agreement of 0.49, 0.47 and 0.46, respectively. From table 6, one can observe a similar (but lower) agreement between the parameters RI, R_{1_d} and R_{1_2} and the reference parameters, with the highest level of agreement being found between RI and SBP (ρ_p = 0.52) and between RI and TPRI (ρ_p = 0.50). The parameter R_{1_2} presented low correlation (${{\rho}_{\text{p}}}\epsilon \left[0.23;0.27\right]$) with all the reference parameters.

The level of agreement between the surrogates presents the best results (i.e. SI) and the SBP are illustrated by the regression plots present in figures 12(a) and (b), for two volunteers. A linear model SBP = A * SI + B was applied on the data extracted from volunteer 7 (figure 12(a)) and volunteer 36 (figure 12(b)). For volunteer 7, the best linear fit corresponds to a negative slope A =  −1.8, which is very far from the unit. For volunteer 36, the best linear fit corresponds to similar slope A =  −2.28. From figure 12 it is possible to observe a strong linear and negative correlation between SI and SBP, especially for lower values of SBP (high values of SI). The higher dispersion SI/SBP points for high SBPs, particularly above 130 mmHg, suggest a higher uncertainty in the SI measurement, resulting from the increasingly overlapping of the forward and reflected waves.

Zoom In Zoom Out Reset image size

The results presented in this section suggest that the extracted parameters can be associated with changes in blood pressure and vascular resistance. The presented results regarding the SI and T_{1_d} parameters for SBP, DBP and MBP are consistent with those found by Millasseau et al (2003) (ρ_p = 0.48/0.68/0.51 for SBP/DBP/MDP), by Padilla et al (2009) (ρ_p = 0.41/0.40/0.44 for SBP/DBP/MDP) and by Brillante et al (2008) (ρ_p = 0.42/0.52/0.51 for SBP/DBP/MDP). Contrarily, the low agreement between RI and blood pressure measurements reported by in (Millasseau et al 2003, Brillante et al 2008, Padilla et al 2009) (${{\rho}_{\text{p}}}$  <  0.18 for SBP, ${{\rho}_{\text{p}}}$  <  0.28 for DBP and ${{\rho}_{\text{p}}}$  <  0.24 for MBP) was not found in our analysis, where the correlation found with SBP, DBP and MBP was above 0.38 for the global dataset and above 0.46 for the tilt positive dataset. The low correlation found between the R_{1_2} and SBP was not consistent with the association between reported by Baruch et al (2011). This discrepancy is mainly explained by the normalization of the PPG pulse performed in the present study, and by the modalities used in the study reported in (Baruch et al 2011). While Baruch et al (2011) used the absolutes amplitudes of the main forward and first reflection waves (assessed from the arterial blood pressure waveform) to estimate R_{1_2}, in the present study the R_{1_2} parameter was assessed from the normalized pulse (assessed from the photoplethysmogram). The lower correlation between SI, T_{1_d} and PP (${{\rho}_{\text{p}}}\approx 0.33$) found in the present study is in agreement with that found by Padilla et al (2009) (${{\rho}_{\text{p}}}=0.24$), but not with the association outlined in (Baruch et al 2011).

The increase in the agreement between the SI, RI and the reference parameters, from the global dataset to the tilt positive dataset, show the importance of these parameters in the assessment of blood pressure and vascular tone changes, especially in cases where acute vasodilatation is present, which is the case of the patients suffering from syncope (tilt positive group).