Heart Rate Variability Assessment Using Time–Frequency Analysis in Hypotensive and Non-Hypotensive Patients in Hemodialysis

Featured Application

In this work, a spectrotemporal evaluation of heart rate variability during hemodialysis shows differences between hypotensive and non-hypotensive patients in the spectral bands related to the sympathetic and parasympathetic nervous system. This methodology could be implemented in a dedicated electronic system for the early detection of intradialytic hypotension events.

Abstract

Intradialytic hypotension occurs in 10–30% of hemodialysis (HD) sessions. This phenomenon affects the cardiovascular system’s functions, which are reflected in the activity of the autonomic nervous system (ANS). To indirectly assess the ANS during HD, we analyzed the mean R–R intervals and the spectral power of heart rate variability (HRV) from 20 end-stage renal disease patients divided into hypotensive and non-hypotensive groups. The spectrotemporal analysis was accomplished using short-time Fourier transform with 10 min epochs of HRV overlapping by 40%. The spectral power was divided into three segments according to high frequency, low frequency, and very low frequency bandwidths and averaged to fit quadratic regression models. The analysis of the mean R–R intervals showed significant differences between the groups (p = 0.029). The power variation over time was significant in each spectral band (p ≪ 0.05). The average power, maximum power, and time when the peak was reached differed for each band and between groups, showing the ability to correctly identify the decompensation of the ANS and discriminate between hypotensive and non-hypotensive patients. Additionally, the changes in the sympathovagal ratio were not significant and very scattered for the hypotensive group (p = 0.23) compared to the non-hypotensive group, where the changes were significant (p ≪ 0.05) and much less scattered.

1. Introduction

The patient population with end-stage renal disease (ESRD) is increasing worldwide [1]. Hemodialysis (HD) is one of the three renal replacement therapies for ESRD laterally with peritoneal dialysis and kidney transplant.

HD is the second most convenient treatment after kidney transplant, but it has severe negative outcomes, as well as high costs and quality assurance (QA) difficulties [1]. The intradialytic hypotension rate is the principal HD difficulty since it is associated with high lethality and a loss of QA. In Mexico, more than 65% of clinics have poor QA figures and report hypotension rates above 30%, leading to unacceptable lethality rates as high as 21% in less than 5 years. This situation has concerned the Mexican Health Care Delivery System since 2011, which has devoted research efforts to find a new method for prognosis to prevent hypotension events in patients undergoing HD [2,3].

The problem is that hypotension events frequently occur without any alert or previous notice since they depend on the cardiovascular (CV) response of each patient when they gradually lose blood fluids and solutes through HD treatment, despite the use of well-established medical standards, such as appropriate ultrafiltration rates and other HD therapy factors [2,4].

According to the Kidney Disease Outcomes Quality Initiative (KDOQI) Guidelines 2015, hypotension is defined as a decrease in the systolic blood pressure by 20 mmHg or of the mean blood pressure by 10 mmHg, accompanied by symptoms such as nausea, vomiting, muscle cramps, and dizziness or fainting [4]. A specific difficulty well established in the literature is that the continuous monitoring of heart rate and mean blood pressure in clinical practice is not sufficient to prognose or prevent the CV system response [5,6,7].

The analysis of heart rate variability (HRV) has already proven its ability to indirectly assess the activity of the autonomic nervous system (ANS). HRV estimates the cyclic modulatory control of the CV system and as a consequence, the link to the blood pressure mechanism. The sympathetic and parasympathetic activity of the ANS is classically measured in the time or frequency domain by using linear tools [8,9]. Specifically, the power spectral analysis of HRV requires a tachogram, which is built using the R–R intervals from a standard electrocardiogram (ECG). The standard is that the total power spectrum from 0 to 0.5 Hz is divided into four bandwidths, identified as ultra-low frequency (ULF, <0.003 Hz), very low frequency (VLF, 0.003–0.04 Hz), low frequency (LF, 0.04–0.15 Hz), and high frequency (HF, 0.15–0.4 Hz) bands. The normal assumption is that each bandwidth is related to different neural controls in the CV system [10]. Thus, the accepted premises are (a) that the HF band is associated with the activity of the parasympathetic nervous system (PNS), (b) that the LF band reflects the activity of the sympathetic nervous system (SNS) [11], (c) that the sympathovagal ratio (LF/HF) measures the proportion of the SNS and the PNS that interacts to control the CV system [12,13], and (d) that the VLF and ULF bands are associated with slow-acting physiological processes, such as circadian rhythm, metabolism, thermoregulation, vasomotor tone, and the renin–angiotensin–aldosterone system (RAAS), which in turn participates in controlling blood pressure and has been closely linked to hypotension events [14,15,16,17].

The spectral power of HRV can be estimated through several mathematical tools [18,19,20], including the fast Fourier transform, autoregressive models, and approximate entropy [21,22]. However, these tools do not provide information related to power variation over time, which is useful in some studies. Therefore, it is important to run a spectrotemporal or time–frequency analysis of HRV. Tools like short-time Fourier transform (STFT), wavelet transform (WT), and the less known Hilbert–Huang Transform (HHT) have been explored for this purpose [23].

The spectral analysis of HRV has been used to identify negative events, such as hypotension and bradycardia, in high-risk cardiovascular patients [24]. It has also been applied to ESRD patients before and after they undergo HD therapy. The results of these works showed higher SNS activity compared to PNS activity at the end of each HD session [25,26,27,28,29,30]. Similarly, Park et al. performed a spectral analysis of HRV over segments acquired at the beginning, middle, and end of HD therapy [31]. The analysis included some time-domain parameters, such as the mean R–R, the standard deviation of all normal to normal (NN) intervals (SDNN), and the root mean square of successive NN differences (RMSSD), as well as blood tests and clinical parameters. Even though including a third point in the spectral analysis improved the assessment of the SNS and PNS during therapy, the approach did not give information on the gradual adaptation of the ANS to the hypovolemic state induced by HD. Importantly, Park’s work showed that HRV can be useful in predicting hypotensive events. In a study published by Huang et al., the authors searched for the association between HRV and the presence of major adverse events related to mortality in patients receiving HD therapy [32]. Their results showed that a decrease in VLF power is related to an increased risk of major adverse events. This work shows that the evaluation of cardiac function through a spectral analysis of HRV has fundamental prognostic importance for these patients.

The hypothesis in this work is that the SNS decreases over time during HD therapy in patients experiencing hypotension events, and that patients who do not show any adverse signs maintain or increase their SNS activity over the duration of HD therapy. To test this hypothesis, we estimated the spectrotemporal power variations of the VLF, LF and HF bands during HD therapy in hypotensive and non-hypotensive patients. We applied the STFT with an epoch twice the size of the standard recommended by the Task Force for HRV Analysis for short-term recordings [13]. The changes in the LF and HF power are related to the adaptive response of the SNS and PNS induced by HD therapy in the CV system [8]. These variations could serve as markers to identify patients in risk of hypotension during HD [31]. Considering the dynamic characteristics of the HRV, it becomes necessary to employ time–frequency analysis tools to track the changes in the spectral power over the intradialysis time. In this work, we employed the STFT because it is a suitable simple-to-use tool that allows to locate the changes in the spectral power over time [33,34].

Unlike other jobs, we included the VLF band in our analysis, since it is considered that the VLF power reflects the activity of the RAAS [16], which participates in the long-term regulation of blood pressure through angiotensin, whose level varies during HD [30].

Previous studies have shown that using 10 min epochs improves the estimation of the VLF spectral power [21,35]. Here, we applied the STFT using this modification in the analysis epoch to study the changes in the power spectral distribution along the intradialysis time, comparing hypotensive and non-hypotensive patients. Along with the spectrotemporal analysis, we also analyzed the mean R–R intervals.

2. Materials and Methods

2.1. Patients and Therapy

A total of 20 ESRD patients participated in this study: twelve female patients with a mean age of 37 ± 14.2 years old and a mean body mass index (BMI) of 21 ± 8.7 Kg/m², and eight male patients with a mean age of 35 ± 14.2 years old and a mean BMI of 22 ± 3.1 Kg/m². All participants received hemodiafiltration (HDF) therapy three times per week at the Hemodialysis Unit at the Instituto Nacional de Cardiología Ignacio Chávez (INCICh) for at least six months before the start of this protocol.

None of the patients were prescribed anti-hypertensive medication, and all of them had stable systolic and diastolic blood pressure (134 ± 22.6 mmHg and 76.7 ± 15.2 mmHg, respectively) at the day of their participation in the study. Comorbidities present in some of the patients were diabetes mellitus (N = 2), arterial hypertension (N = 4), and hepatitis C (N = 4). All patients were volunteers and signed their informed consent before participation. The protocol was approved by the Institutional Ethical Committee of the INCICh (number 17-1014) and was conducted in accordance with the Declaration of Helsinki.

The HDF sessions lasted at least 3 h and were performed in the first morning shift using the same Fresenius 4008 hemodialysis machine (Fresenius Medical Care). The ultrafiltration rate was standardized to 800 mL/h for a total volume extraction of 2400 mL. Patients were asked to follow dietary recommendations to try not to exceed this interdialytic volume gain between two consecutive sessions. Additional measures were taken in case the patient gained a different interdialytic volume than the target volume. If a larger volume gain happened, then the HDF session was extended to reach the pre-established dry weight, but the ECG was acquired only over the first three hours. If the volume extraction was short, then an additional bolus was injected during the HDF session to reach the 2400 mL established for the study. Hemodialysis therapy was carried out in a sitting position. From the 20 patients, 8 experienced an intradialytic hypotension episode during their therapy (7 women and 1 man). In these cases, the Trendelenburg maneuver was applied in combination with a decrease in the ultrafiltration rate to stabilize the patient.

2.2. ECG Recording System

A 7-lead Holter device (model DMS300-7, DM Systems Co., Ltd., Changping, Beijing, China) with a sampling rate of 128 samples/second was used to record the ECG during HDF. The leads were placed in an orthogonal configuration (Frank leads) to acquire three channels simultaneously, and the data acquisition began at least 10 min before the HDF therapy started. At the end of the therapy, the recordings were downloaded using CardioScan^® software for storage in a desktop computer and further offline processing and analysis.

2.3. ECG Processing and Generation of the HRV Signal

According to the report yielded by CardioScan, all ECGs had less than 1% premature beats. The recordings were resampled to 256 samples/s using linear interpolation with MATLAB R2016 (MathWorks Inc., Natick, MA, USA) before filtering by moving averages and generating the HRV recording. To upsample the ECG-signal before R-wave detection is highly advisable [13], and a minimum of 250 samples/s is recommended to not compromise the high frequency content [36]. The three acquired channels were visually inspected to choose the one with less distortion for further processing. The R-wave detection, visual inspection and manual correction of misidentified peaks, the generation of the R–R tachograms, and the regular resampling of the tachograms to 4 samples/s using cubic spline interpolation were done with the Kubios HRV Analysis Software [37]. Regular sampling is mandatory in STFT given the irregular location of the detected R-waves. As an example of the generation of the HRV signal, Figure 1 shows a 10 s segment of the filtered ECG signal and the trace corresponding to the QRS complex detection (panel a), the identification and location of the corresponding R-waves (panel b), and the irregular R–R intervals obtained along altogether with the R–R tachogram resampled at 4 samples/s (panel c).

2.4. Short-Time Fourier Transform Theory

Fourier transform (FT) is a mathematical tool used to perform the linear decomposition of a signal into its constituent spectral components. The resulting decomposition provides spectral information of the analyzed signal but no time localization of the spectral content [18,19]. To apply FT, it is necessary the signal to be stationary.

Short-time Fourier transform is a Fourier-related transform able to deal with the non-stationarity characteristics of a signal, albeit in a limited way due to the length of the signal segment to be analyzed. The results yield information about the time localization of the spectral components, which are computed by performing a Fourier Transform over overlapped windowed segments of a signal to reduce spectral leakage [33]. The continuous form of the STFT is expressed in Equation (1):

$$STFT\left\{x\left(t\right)\right\}=X\left(\tau ,f\right)={{\displaystyle \int }}_{-\infty }^{+\infty }x\left(t\right)w\left(t-\tau \right)e^{-j2\pi ft}dt$$

(1)

where $x\left(t\right)$ is a time-varying signal, $w\left(t-\tau \right)$ is a window function shifted by the parameter $\tau$, $e^{-j2\pi ft}$ is the Euler identity, and $X\left(\tau ,f\right)$ is the Fourier transform of the windowed segment of $x\left(t\right)$. The discrete form of the STFT is expressed in Equation (2):

$$STFT\left\{x\left(n\right)\right\}=X\left(m,k\right)={\displaystyle \sum }_{n=0}^{L-1}x\left(n\right)w\left(n-m\right)e^{-j2\pi n\frac{k}{L}}$$

(2)

where $x\left(n\right)$, $w\left(n-m\right)$, $e^{-j2\pi n\frac{k}{L}}$, and $X\left(m,k\right)$ represent the same variables as those in the continuous form but dependent on the discrete variable $n$ and $k$.

The results of the transformation can be displayed using a spectrogram, which provides a graphical representation of the spectral power distribution over time. This is obtained by squaring the absolute value of Equation (2):

$$S_x^w\left(m,k\right)={\left|X\left(m,k\right)\right|}^2={\left|{\displaystyle \sum }_{n=0}^{L-1}x\left(n\right)w\left(n-m\right)e^{-j2\pi n\frac{k}{L}}\right|}^2$$

(3)

The spectrogram obtained with the STFT is distorted by the spectral leakage resulting from the windowing operation and the Heisenberg’s uncertainty principle. In signal processing, this principle states that an accurate estimation of both time and frequency cannot be obtained simultaneously. This implies that a wide window gives better frequency resolution but poor time resolution, while a narrower window gives good time resolution but poor frequency resolution.

Since the length of the window is fixed and therefore, cannot be modified during the transformation, a balance must be determined a priori between the desired time and frequency resolutions according to the needs of the analysis. To determine the frequency resolution, we took into account the lower frequency limit of the VLF band, which is 3 mHz. To represent one full cycle of this frequency, we need 333.33 s. However, for nondeterministic signals (as is the case of the HRV), it is recommended to use analysis epochs longer than the slowest spectral component [13]. In a previous work, it is shown that the spectral power estimation of a multicomponent signal decreases as the epoch increases with respect to the slowest component [21]. Based on that work, we decided to use 10 min epochs. This decision also considered that the number of segments to be analyzed during HD decreases as the epoch increases, and becomes less representative of the dynamics of the estimated spectral power. Once we determined the time resolution, we calculated the frequency resolution obtained based on the uncertainty principle, which states that:

$$\Delta t\Delta f\ge \frac{1}{2}$$

(4)

where $\Delta t$ is the time resolution and $\Delta f$ is the frequency resolution. For $\Delta t=600s$, $\Delta f$ is at least 0.833 mHz, which is lower than the lower limit of the VLF band.

2.5. HRV Processing

Before applying STFT, the R–R tachograms were extended at both ends by mirroring the first and last 5 min. This preprocessing is mandatory because the FT is centered in the sliding window. If this signal extension was not made, we would lose the information present at the beginning and the end of the tachograms. To diminish spectral leakage, the Hanning window was used in the transformation [38].

To obtain the maximum number of independent analysis points, each 3 h R–R tachogram was split into 30 segments of 10 min with a 40% (4 min) overlap between adjacent segments. This percentage was determined by running an autocorrelation analysis each time the overlap increased by 10%, from 0% to 100%. Independence was found with a 0% to 40% overlap between consecutive segments. On the contrary, no independence was found with a 50% to 100% overlap.

2.6. Statistical Analysis

The correlation of the mean spectral power with time was assessed using quadratic regression models. The regression coefficient (R), the determination coefficient (R²), and the p-value were also calculated from each model for each spectral band. Statistical significance was evaluated using an F-test.

Along with the spectral analysis, the R–R intervals were also analyzed. Two series were obtained by averaging the R–R intervals in each of the 30 segments and for each group of hypotensive and non-hypotensive patients. A t-test was run on these two series to determine if there was a significant difference between them. The normal distribution of the R–R intervals was tested using the Shapiro–Wilk normality test. Correlation with time was also assessed using quadratic regression models. All results were considered statistically significant for p ≤ 0.05. The statistical analysis was done using Statistic^® ver. 8 (TIBCO Software Inc., Palo Alto, CA, USA).

We verified the achieved statistical power of the quadratic regression analyses with the program G*Power (version 3.1.9.7). For each regression model, we used the obtained determination coefficient (R²) to calculate the size effect, and then we estimated the achieved statistical power for the model with two predictors, a sample size N = 30, and an error probability α = 0.05.

3. Results

According to the results obtained with the Shapiro–Wilk test, the two mean R–R series are normally distributed (p > 0.05). Figure 2 shows the mean R–R intervals plotted as a function of time for each group of patients. The first mean value obtained from the R-R intervals is placed at 10 min, and the last mean value is at 184 min. Adjacent points are separated by 6 min according to the maximum overlap described in the HRV processing section to maximize the number of analysis points. The mean and standard deviation obtained were 814.38 ± 20.40 ms for the group of hypotensive patients and 829.23 ± 30.05 ms for the group of non-hypotensive patients. The t-test showed significant differences between the groups (p = 0.029). The regression models obtained were y = 838.36 − 0.24t − 0.00005t² (ms) (p ≪ 0.05, R = 0.65, and R² = 0.42) for the hypotensive group and y = 757.56 + 2.18t − 0.012t² (ms) (p ≪ 0.05, R = 0.95, and R² = 0.90) for the non-hypotensive group. The lowercase t in the models represents time.

From the time–frequency analysis, the spectrograms obtained show clear differences in the power distribution for the hypotensive group compared to the non-hypotensive group. Figure 3 shows the spectrograms obtained from two patients—one from each group—for the visualization purposes. Although the differences are spread over the time–frequency plane, a greater contrast was observed in the HF band, where the power was smaller for the hypotensive patients compared to the non-hypotensive patients.

Figure 4 shows the mean spectral power obtained from the VLF, LF, and HF bands for the two studied groups of patients, along with their corresponding regression models and confidence intervals. The VLF power is lower for the hypotensive group compared to the non-hypotensive group (Figure 4, panel a). The obtained regression models were y = 557.04 + 10.6t − 0.050t² (ms²/Hz) (p ≪ 0.05, R = 0.81, and R² = 0.66) and y = 514.88 + 20.87t − 0.092t² (ms²/Hz) (p = 0.001, R = 0.53, and R² = 0.28), respectively. The models reached their maximum at 107 min (hypotensive group) and 113 min (non-hypotensive group).

For the LF power, the models obtained were y = 475.13 + 7.88t − 0.027t² (ms²/Hz) (p ≪ 0.05, R = 0.77, and R² = 0.59) for the group of non-hypotensive patients and y = 556.29 + 7.27t − 0.052t² (ms²/Hz) (p ≪ 0.05, R = 0.68, and R² = 0.46) for the group of hypotensive patients (Figure 4, panel b). According to these models, the maximum is reached at 145 min for the first group and at 70 min for the second group. This causes the LF power to decrease first in the second group and reach a smaller value at the end of the HD therapy compared to the first group.

The results for the HF power and the two groups of patients are shown in Figure 4 (panel c). The model corresponding to the non-hypotensive group is y = 153.67 + 6.76t − 0.030t² (ms²/Hz) (p ≪ 0.05, R = 0.87, and R² = 0.76), and the one corresponding to the hypotensive group is y = 157.17 + 0.84t − 0.009t² (ms²/Hz) (p ≪ 0.05, R = 0.69, and R² = 0.48). The maximums of the models were reached at 105 min for the first group and at 48 min for the second group.

All times when the maximums were reached were calculated using the corresponding quadratic model.

The results of the sympathovagal ratio for the two groups, as well as the respective regression models and confidence intervals, are plotted in Figure 5. The model for the group of non-hypotensive patients is y = 3.84 − 0.03t + 0.0001t² (p ≪ 0.05, R = 0.69, and R² = 0.48), while that for the group of hypotensive patients is y = 6.33 + 0.04t − 0.0002t² (p ≪ 0.05, R = 0.31, and R² = 0.10).

Table 1. Quadratic regression models, statistical parameters achieved and statistical power for the hypotensive group. Variable t represents time. p < 0.05 was considered significant.

Variable	Quadratic Model	R	R2	p	Achieved Statistical Power
Mean R–R	838.36 − 0.24t − 0.00005t2	0.65	0.42	≪0.05	0.98
VLF	557.04 + 10.60t − 0.050t2	0.53	0.28	0.01	0.83
LF	556.29 + 7.27t − 0.052t2	0.68	0.46	≪0.05	0.99
HF	157.17 + 0.84t − 0.009t2	0.69	0.48	≪0.05	0.99
LF/HF	6.33 + 0.04t − 0.0002t2	0.31	0.10	0.23	0.32

and

Table 2. Quadratic regression models, statistical parameters and achieved statistical power for the non-hypotensive group. Variable t represents time. p < 0.05 was considered significant.

Variable	Quadratic Model	R	R2	p	Achieved Statistical Power
Mean R–R	757.56 + 2.18t − 0.011t2	0.95	0.90	≪0.05	1.00
VLF	514.88 + 20.87t − 0.092t2	0.81	0.66	≪0.05	0.99
LF	475.13 + 7.88t − 0.027t2	0.77	0.59	≪0.05	0.99
HF	153.67 + 6.76t − 0.030t2	0.87	0.76	≪0.05	1.00
LF/HF	3.84 − 0.03t + 0.0001t2	0.69	0.48	≪0.05	0.99

summarize the quadratic models, the statistical parameters, and the achieved statistical power for the mean R–R intervals; the VLF, LF, and HF powers; and the LF/HF ratio. All models except one achieved a statistical power above 0.80.

4. Discussion

The objective of this work was to indirectly assess the modulatory activity of the ANS in ESRD patients through a spectrotemporal analysis of HRV and to evaluate the behavior of the mean R–R intervals (both during hemodiafiltration). The participants included patients who experienced an intradialytic hypotension event and patients who did not. This characteristic allowed us to split the participants into two groups and study the differences in the modulatory response of the SNS and PNS when this adverse event occurs. Furthermore, we also showed the feasibility of using 10 min epochs overlapping by 40% to perform the analyses and reduce the errors in the spectral estimation of VLF power.

The results obtained with the mean R–R intervals show visual and statistical differences between the groups, which indicate that it is possible to use this measure to discriminate between hypotensive and non-hypotensive patients.

To assess the spectrotemporal distribution of power, we used the STFT due to its ability to operate over short segments of a signal that can be considered stationary and linear. This characteristic diminishes the difficulties found in analyzing HRV because of the dynamics present in the CV system. Other advantages of the STFT are its low computational load compared to other time–frequency analysis tools [23,33] and its lower amount of spectral leakage resulting from the use of a single epoch to perform the analysis [38], which allows the identification of large spectral components present in noisy signals.

The time–frequency analysis was performed using 10 min segments of HRV. According to the guidelines published in 1996 by the joint American and European Task Force for HRV Analysis [13], this spectral analysis can be performed in either short-term recordings 2 to 5 min in length or long-term recordings of 24 h in length. When short-term recordings were used, the guidelines indicated that it is possible to estimate the power spectrum in the LF, HF, and VLF bands; however, it is also advised in the guidelines that the results obtained for the VLF power are not reliable; hence, its estimation must be avoided. For long-term recordings, the guidelines were not restricted to estimating the power in any of the three spectral bands. Along with the standardization of the length of the analyzed segments, the Task Force also recommends removing the mean and trend of the signal before performing a spectral analysis to minimize the spectral contribution below the VLF band to 0 Hz.

Since it is not reliable to estimate the VLF power in recordings ≤5 min, Becerra et al. [21] evaluated the errors obtained when estimating the spectral power distribution of a signal when the length of the analyzed segment varies from 300 s to 4000 s. The results show that the error decreases as the length increases. According to the publication, to decrease the errors in the VLF power estimation to slightly less than 20%, the length of the segment must be increased to 10 min. Following these results, we performed the STFT over 10 min segments.

The results of the spectrotemporal analysis show that all powers behave similarly in the three bands. We observed that any spectral power can be modeled by a vertical parabola opened downward (Figure 4). It can also be seen that the regression models of hypotensive patients had higher mean values than those of non-hypotensive patients and that the time when the maximum power was reached (obtained with the respective quadratic models) is different for each band and depends on the presence or absence of a hypotension event. Moreover, the plots of each band show that the models for both groups have very close starting power values. The statistics of all fitted models yielded significant p-values for their correlation with time.

Despite the high similarity between the two quadratic models for VLF power in the independent term and moment when the maximum is reached, the main difference between them is the rate of change, which is smaller for the hypotensive group than for the non-hypotensive group. This caused the maximum power reached by the first group to be smaller than the maximum power of the second group, although there was only 6 min of difference between them. Slow changes in the VLF power might suggest the decompensation of the CV system, which is strongly linked to an increased risk of developing major adverse cardiovascular events [32]. Although the LF and HF powers received more attention in the spectral analysis, it is important to follow changes in the VLF band because such changes are associated with the RAAS [16], which has a strong bidirectional interaction with the autonomic nervous system [39].

For the LF power, the maximum was reached during the first half of HD therapy (at 70 min) for the group of hypotensive patients. Here, the activity of the sympathetic system responded normally by increasing the heart rate. However, for the remaining time of the HD therapy, the LF power decreased slowly, while the heart rate continued to increase, showing behavior contrary to what was expected. For the group of non-hypotensive patients, this value reached its maximum at 145 min, and during the remaining 35 min of therapy, the LF power decreased minimally, favoring a lower reduction in heart rate and sustained arterial pressure for a longer time.

The HF power was larger for the group of non-hypotensive patients, and the maximum value of the model occurred almost during the third hour of therapy (116 min), while for the group of hypotensive patients, the maximum value was reached during the first hour of therapy. Further, the HF power of the hypotensive patients increased very little compared to that of the non-hypotensive patients, which suggests continuous attenuation in the vagal modulatory mechanism during the therapy.

All the quadratic models obtained for the spectral powers yielded R > 0.52 and R² > 0.27, which indicates that the model is correctly able to describe the behavior of HRV during HDF, with statistical power greater than 0.8 in all models except one.

The results we obtained are similar to those published by Park et al. [31], who reported an analysis of R–R intervals in short-term recordings acquired at the beginning, middle, and end of HD therapy. Even though we used more analysis points obtained over longer epochs (30 points obtained from overlapping segments of 10 min each versus three points from non-overlapping segments of 5 min segments each), we found significant differences (p < 0.05) between hypotensive and non-hypotensive patients for the mean R–R intervals obtained both in the middle and at the end of the therapy. The results of the spectral analysis agree in terms of behavior and significant differences for each of the three bands evaluated. An important characteristic to mention is that both our studies used absolute spectral power (ms²/Hz) to compare the two groups of patients.

On the contrary, the results of the spectral analysis published by Barnas et al. [8] show that when a hypotension event occurs, the LF power decreases, and the HF power increases (both greatly) compared to their respective values in the middle of therapy; for those patients who did not experience a hypotension event, the LF power increases, and the HF decreases at the end of the therapy compared to their respective values in the middle of therapy. This difference with our results might be attributable to the normalization procedure the authors used before comparing the spectral power from hypotensive and non-hypotensive patients. Both studies also had a similar number of patients (20 in this work and 19 in Barnas’ research).

Another interesting result we found is that the LF/HF ratio is very scattered in hypotensive patients, which explains the low correlation and determination coefficients obtained with this model. The high dispersion of these measurements could be associated with dysfunction in the autonomic regulation of these patients, which greatly attenuates HF’s power. On the contrary, the LF/HF ratio shows much less dispersion in non-hypotensive patients and can be well adjusted to a quadratic regression model.

It is worth noting that each group of patients has different average LF/HF ratios during therapy. For hypotensive patients, this ratio is close to 8, while for non-hypotensive patients, it is close to 3.

The results presented in this work show that the decompensation of the CV system during HD therapy, leading to a hypotension event in some patients, is reflected by the modulatory activity of the SNS and PNS through changes in the magnitude of the spectral power in the VLF, LF, and HF bands over time. These differences can be evaluated through the application of the STFT, by allowing to know the dynamics of the HRV through time.

Two points that stand out in the time–frequency analysis with the STFT are related to the length of the window chosen for performing STFT since this length determines the time–frequency resolution, and the window function was employed because the spectral leakage depends on this function. In our application, the length of the window was determined by the size of the epoch required to decrease the error in the VLF power estimation to at least 20%, and the window function chosen was the Hanning or Raised–Cosine window due to its minimal spectral leakage [38].

A limitation of our study is the relatively small size of our population (N = 20, 8 hypotensive versus 12 non-hypotensive). Despite this small number of participants, the differences found in the spectral analysis are clear enough to discriminate patients who experienced hypotension events from those who did not.

With regard to the time–frequency analysis tools, we acknowledge that there are other alternatives to the STFT, including the wavelet transform, the Hilbert–Huang transform [23], and the synchrosqueezed transform [40]. Further studies are required to compare the estimation of the spectrotemporal power of HRV during HD with the other time–frequency analysis tools, and to assess their effectiveness. In this work, we chose the STFT because of its simplicity and low computational cost. These characteristics make the STFT a good candidate for further implementation in a dedicated electronic system, to run the analysis online during hemodialysis therapy using the type of model presented in this research.

5. Conclusions

The HRV signal (R–R tachogram) is a time series that is easily obtained from an electrocardiogram. The mean R-R interval calculated using 10 min epochs is useful for discriminating between hypotensive and non-hypotensive patients.

The spectrotemporal analysis of HRV is a more laborious method to discriminate between hypotensive and non-hypotensive patients but offers more information related to the autonomic nervous system over time. The spectrograms obtained with STFT show differences between these two groups of patients. By averaging the spectral power obtained from each band and group of patients, we can build quadratic regression models to evaluate the correlation with time.

The power variation during the therapy has a similar shape for the three bands and the two groups of patients, albeit important differences can be seen. For non-hypotensive patients, the power in the three bands was larger, and reached the maximum and decreased during the second half of the therapy. Even though the decreasing the corresponding power in the three bands were larger at the end of the HD compared to the beginning. The LF power in this group kept increasing almost during the whole therapy. This behavior along with the behavior of the HF band suggests a sustained activity of the SNS and PNS that allows the regulation of the CV system and avoids the advent of a hypotension event.

For the group of hypotensive patients, the maximum in each band was reached at a different moment of the HD therapy. Special attention should be placed on the LF and HF powers, whose values were smaller at the end compared to the beginning. These findings suggest an imbalance of the ANS, that leads to a hypotension event, in agreement with our hypothesis.