Monitoring of Serum Potassium and Calcium Levels in End-Stage Renal Disease Patients by ECG Depolarization Morphology Analysis

Objective: Non-invasive estimation of serum potassium, [K+], and calcium, [Ca2+], can help to prevent life-threatening ventricular arrhythmias in patients with advanced renal disease, but current methods for estimation of electrolyte levels have limitations. We aimed to develop new markers based on the morphology of the QRS complex of the electrocardiogram (ECG). Methods: ECG recordings from 29 patients undergoing hemodialysis (HD) were processed. Mean warped QRS complexes were computed in two-minute windows at the start of an HD session, at the end of each HD hour and 48 h after it. We quantified QRS width, amplitude and the proposed QRS morphology-based markers that were computed by warping techniques. Reference [K+] and [Ca2+] were determined from blood samples acquired at the time points where the markers were estimated. Linear regression models were used to estimate electrolyte levels from the QRS markers individually and in combination with T wave morphology markers. Leave-one-out cross-validation was used to assess the performance of the estimators. Results: All markers, except for QRS width, strongly correlated with [K+] (median Pearson correlation coefficients, r, ranging from 0.81 to 0.87) and with [Ca2+] (r ranging from 0.61 to 0.76). QRS morphology markers showed very low sensitivity to heart rate (HR). Actual and estimated serum electrolyte levels differed, on average, by less than 0.035 mM (relative error of 0.018) for [K+] and 0.010 mM (relative error of 0.004) for [Ca2+] when patient-specific multivariable estimators combining QRS and T wave markers were used. Conclusion: QRS morphological markers allow non-invasive estimation of [K+] and [Ca2+] with low sensitivity to HR. The estimation performance is improved when multivariable models, including T wave markers, are considered. Significance: Markers based on the QRS complex of the ECG could contribute to non-invasive monitoring of serum electrolyte levels and arrhythmia risk prediction in patients with renal disease.


Introduction
Patients with chronic kidney disease (CKD), particularly at the most severe stages, are at high risk for life-threatening arrhythmias and sudden cardiac death [1][2][3][4], and millions of them die every year due to lack of access to affordable treatment. Between five and seven million end-stage renal disease (ESRD) patients need renal replacement therapy worldwide [5]. ESRD patients undergoing hemodialysis (HD) frequently have serum electrolyte levels outside normal ranges, which can increase the mortality risk [1][2][3][4].
Non-invasive monitoring of electrolyte levels can be useful for risk prediction and the triggering of early warnings. Electrocardiogram (ECG) ventricular depolarization and repolarization are, in particular, altered by variations in serum potassium ([K + ]) and calcium ([Ca 2+ ]) [1,2,[6][7][8]. ECG markers have been proposed for the continuous monitoring of these two electrolyte levels, which could facilitate timely therapies for ESRD patients.
Most ECG-based estimators of [K + ] and [Ca 2+ ], such as those based on the QT interval or on T wave amplitude, slope and slope-to-amplitude ratio, are based exclusively on ventricular repolarization and commonly refer to a specific time interval, to an amplitude at a given time point or to a particular portion of the T wave [9][10][11][12][13][14][15]. Although the QT interval has long been used to monitor CKD patients during HD [16][17][18][19][20][21], contradictory findings have been presented, with many studies reporting QT prolongation [18,19,[21][22][23] and others showing QT shortening or no effects during HD [24,25]. Other recent studies have investigated changes in markers of sympathetic activity-related T wave instability during HD [26] and in model-based descriptors of T wave morphology in the interdialytic interval [27]. These studies, however, have not established a tight correlation between the proposed markers and [K + ], thus limiting their possibilities for ambulatory [K + ] monitoring. In previous studies, we used nonlinear dynamics and time-warping techniques to characterize changes in the whole T wave at varying [K + ] and [Ca 2+ ] in patients and simulated ECGs [28][29][30][31][32][33]. We found a strong relationship between [K + ] and T wave linear and nonlinear features.
Assessment of ventricular depolarization for serum electrolyte estimation has been less explored. Similarly to QT studies, research on QRS complex duration has rendered inconsistent results, with some works reporting widened QRS complexes [1,34] and others reporting narrowed QRS complexes at high [K + ] [2,35]. Other works have assessed the time voltage area, amplitude and sine wave shape of the QRS complex, but limitations in terms of their significance or their dependence on blood volume have been acknowledged [14,[36][37][38]. None of these studies have, however, quantitatively characterized overall variations in the morphology of the QRS complex during HD. We hypothesized that QRS morphology could provide complementary information on [K + ] and [Ca 2+ ], in addition to that provided by repolarization.
The aim of this study was to quantify changes in the QRS amplitude, duration and morphology, the latter both in the time and amplitude domains, at varying [K + ], [Ca 2+ ] and heart rate (HR) in ESRD patients. Univariable and multivariable electrolyte estimators, including novel QRS morphological markers in combination with already proposed T wave markers [28][29][30]39], were devised, and the contribution of depolarization analysis to [K + ] and [Ca 2+ ] monitoring was assessed.

Materials
Forty-eight-hour 12-lead ECGs with 3.75 µV resolution and a sampling frequency of 1 kHz (H12+, Mortara Instruments, Milwaukee, WI, USA) were acquired from 29 ESRD patients of Hospital Clínico Universitario de Zaragoza (HCUZ). ECG acquisition started 5 minutes before HD and lasted for 48 hours (Figure 1, bottom blue line), with patients in the supine position. Concurrently, six blood samples were taken during and after the HD session: the first one at HD onset and the next three samples every hour during the HD session ( Figure 1, h 0 to h 3 in red). The fifth blood sample (h 4 ) was obtained at the end of the HD (minute 215 or 245, depending on the patient) while the sixth blood sample was taken after 48 hours (h 48 ), immediately before the next HD session. Serum [K + ] and [Ca 2+ ] were measured at those six time points, as shown in Figure 1 [28][29][30], using a Cobas 6000 c501 analyzer (Roche Diagnostics, Germany) by an indirect ion selective electrode method. We defined another time point for ECG analysis, corresponding to a segment taken two minutes before the end of HD (minute 213 or 243, depending on the patient), which we denoted by h − 4 . The [K + ] value at h − 4 was assumed to be the same as at h 4 , as the time difference between these two segments was just two minutes. All patients signed informed consent. The study protocol was approved by the Research Ethics Committee of Aragón (CEICA, ref. PI18/003, 25 January 2018). Table 1 shows the population characteristics.

Methods
A flow chart showing all the stages of the signal processing starting with 12-lead ECG pre-processing, followed by principal component analysis (PCA) transformation, computation of mean warped QRS (MWQRS) complexes and QRS morphological markers, and finishing with the estimation of [K + ] and [Ca 2+ ] is shown in Figure 2.

ECG Pre-Processing
A band-pass filter (0.5 to 40 Hz) was used to remove baseline wander, muscular noise and powerline interference from measured ECG signals. A wavelet-based single-lead delineation method [40] was used for QRS detection and wave delineation. Spatial principal component (PC) analysis was performed on the complete ECG of the eight independent leads [41] to enhance QRS complex energy. PCs were obtained from the eigenvectors of the 8 × 8 inter-lead auto-correlation matrix of QRS complexes extracted from a stable (in terms of HR) 10-min ECG segment at the end of the HD session to emphasize the QRS complex components. This segment corresponded to the time when the patient was discharged from hospital with restored serum [K + ] and [Ca 2+ ]. The ECG recording from each time point during and after HD was subsequently projected onto the direction of the first PC. The QRS complexes from each time point in the first PC were delineated [40] to mark their onsets, peaks and ends.

Computation of Mean Warped QRS Complexes
MWQRS complexes, which are optimal representative averages in both temporal and amplitude domains, were calculated from the two-minute ECG segment at the end of each HD hour following the methodology described in [39]. The two-minute ECG segment was short enough to maintain the assumption of stability for both electrolyte and HR values [30].
Only QRS complexes in each analyzed two-minute window that presented the dominant polarity were considered for MWQRS calculation, with the polarity defined as: where f (n) represents the QRS complex under analysis. An average of 92% of QRS complexes in each analyzed two-minute ECG segment were found to present such dominant polarity (see Figure S1 in Supplementary Material). An initial MWQRS was computed after aligning the QRS complexes having the dominant polarity with respect to their gravity center so that the calculated MWQRS was not affected by potential outlier QRS complex polarities and morphologies [39]. Outliers, defined as those QRS complexes having a Spearman's correlation coefficient with an initial MWQRS lower than 0.98, were rejected. The final MWQRS was computed from the remaining QRS complexes. The QRS descriptors defined in the next sections were computed from MWQRS complexes for ECG segments during and after HD (h 0 , h 1 , h 2 , h 3 , h − 4 , h 4 and h 48 ).

QRS Duration and Amplitude Markers
For each MWQRS, the following descriptors were computed: • QRS w , which represented the QRS width calculated from QRS onset to end (expressed in ms) [40]. • QRS a , which represented the QRS amplitude calculated from the minimum to maximum amplitude of the QRS complex (expressed in mV).

QRS Morphology Markers
Morphology-based QRS descriptors were computed using the time-warping methodology, as previously described [39]. For each patient, a reference QRS complex, f r (t r ), was calculated from the MWQRS at the end of the HD session, the time when the patient was discharged from hospital with restored serum ion levels. Representative QRS complexes for each HD stage, f s (t s ), were calculated from the MWQRS at two-minute windows at the end of each hour during the HD session and 48 h after it. f s (t s ) is expressed as f s (t s ) = [ f s (t s (1)), ..., f s (t s (N s ))] and the reference QRS complex as f r (t r ) = [ f r (t r (1)), ..., f r (t r (N r ))] . The vectors t r = [t r (1), ..., t r (N r )] and t s = [t s (1), ..., t s (N s )] are the uniformly sampled time vectors corresponding to the QRS complexes f s and f r , respectively. N r and N s represent the total duration, in samples, of t r and t s , respectively. Figure 3a shows an example of f r and f s , with their respective time domains, t r and t s .
Let γ(t r ) be the warping function that relates t r and t s such that f s (γ(t r )) denotes the time-domain warping of f s (t s ) using γ(t r ). The square-root slope function (SRSF) transformation was used to find the optimal warping function [39]. This transformation was defined as: The optimal warping function was determined as the one minimizing the amplitude difference between the SRSF of f r (t r ) and f s (γ(t r )): A dynamic programming algorithm was used to obtain the function γ * (t r ) that optimally warped f r (t r ) into f s (t s ). This function is shown in Figure 3d. The warped QRS complex, f s (γ * (t r )), is shown in Figure 3b, together with the reference QRS complex, f r (t r ).
The descriptor, d u w,Q , was used to quantify the level of warping required to optimally align the QRS complexes f s (t s ) and f r (t r ): The amplitude descriptor, d a,Q , was computed from the area contained between f r (t r ) and f s (γ * (t r )) normalized by the L2-norm of f r (t r ), thus quantifying amplitude differences after time warping the two QRS complexes: where ) was used to account for the sign. The above described marker d u w,Q accounted for both linear and nonlinear warping required to fit the two QRS complexes in the time domain. The nonlinear component marker d NL w,Q was quantified as: where γ * l (t r ) (green line in Figure 3d) was derived by linearly fitting γ * (t r ) using the least absolute residual method.
The nonlinear amplitude marker d NL a,Q was defined by computing the L 2 norm of the difference between L 2 -normalized versions of f r (t r ) and f s (γ * (t r )): Figure 3a shows f r and f s , referring to their respective time domains, t r and t s . The warped QRS complex (red), f s (γ * (t r )), is shown in Figure 3b, together with the reference QRS complex (blue), f r (t r ).
The analyzed warping-based QRS descriptors, computed from MWQRS during and after HD stages, h i , included: • d u w,Q , which represented temporal variations in QRS morphology (expressed in ms), • d a,Q , which represented amplitude variations in QRS morphology (expressed as a percentage), • d NL w,Q , which represented nonlinear temporal variations in QRS morphology (expressed in ms), • d NL a,Q , which represented nonlinear amplitude variations in QRS morphology (expressed as a percentage).

Statistical Analysis
Pearson correlation coefficients (r) were computed to assess the strength of the linear relationship between [K + ], [Ca 2+ ] or RR and each of the investigated QRS descriptors in each patient individually.
A Wilcoxon signed-rank test was applied to test for significant differences (p < 0.05) w,Q and d NL a,Q between consecutive HD stages. The reason for using a non-parametric statistical test was that, according to the Shapiro-Wilk test, the data distributions were not normal.
To test whether the Pearson correlation coefficient between each QRS marker and [K + ], [Ca 2+ ] or RR was significantly different from 0, a Student's t-test was used after transforming the statistical distribution of r into a normal distribution using Fisher's z transform [42]. All statistical analyses were performed using MATLAB version R2020b for Windows (MathWorks Inc., Novi, MI, USA).

Uni-and Multivariable Estimation of [K + ] and [Ca 2+ ]
To estimate [K + ] and [Ca 2+ ], univariable and multivariable estimators were developed. Of all the analyzed QRS descriptors, d u w,Q was selected to build univariable estimators because it presented high median absolute Pearson correlation with [K + ] and [Ca 2+ ] and a relatively low IQR range, particularly for [K + ]. Univariable estimators were additionally built using the repolarization descriptor d u w,T , which had shown strong correlation with electrolyte levels [29,30]. d u w,T was calculated analogously to d u w,Q but for the T wave instead of the QRS complex. Multivariable estimators were tested using both d u w,Q and d u w,T . The univariable [K + ] u ([Ĉa 2+ ] u , respectively) and multivariable [K + ] m ([Ĉa 2+ ] m , respectively) estimators were of the following form: and where The performance of the estimators was tested by using the leave-one-out cross-validation approach.
For the univariable estimators, the coefficient vector T was estimated as [43,44]: with X = j T x T b . Three different types of estimators were considered, namely stagespecific, patient-specific and global estimators, as described in the following paragraphs. The definition of j T , x T b and y was different for each type of estimator: • For an HD stage-specific (S) estimator, which estimates the electrolyte level at stage h i of a given patient q from the marker's values of the remaining patients at that stage, the vectorβ was calculated from the vector j = 1 1 · · · 1 of dimension 1×Q, with Q being the total number of patients minus 1, the vector containing the values of the marker b = d u w,Q or d u w,T at stage i from patients 1, . . . , Q (all except for patient q) and the vector y = [   T was estimated using Equation (11), The relative error, e r , with respect to the range of electrolyte measurements was computed as The mean absolute error was computed by averaging the absolute value of e defined in Equation (12). The root mean square error was computed by taking the root mean square of e defined in Equation (12).
To assess the agreement between actual and estimated [K + ] and [Ca 2+ ] values, Bland-Altman analysis was performed [45], in which the differences between the actual and estimated [K + ] and [Ca 2+ ] were plotted as a function of their averages, for all patients at all HD time points (see Supplementary Material). [K + ] and [Ca 2+ ] estimation accuracy using our investigated markers was also compared with those of the previously proposed markers T S/A [10,11] and T S/ √ A [12] (see Supplementary Material): • T S/A represented the ratio between the maximal downward slope (in absolute value) and the amplitude of the T wave [10,11]. • T S/ √ A represented the ratio between the maximal downward slope (in absolute value) and the square root of the amplitude of the T wave [12].
It should be noted that estimations of [K + ] and [Ca 2+ ] could not be performed at the end of the HD session (h 4 ) since the morphological QRS complex markers were zero by definition due to the fact that the reference was taken at that time stage. Therefore, we defined an extra time point h − 4 just before the HD end (reference) so as to improve the estimation accuracy based on one additional HD point.

Characterization of QRS Complex Changes during and after HD
The top panels in Figure 4 Figure S2 in the Supplementary Material shows the relative errors, e r . Table 3 shows actual and estimated [K + ] and [Ca 2+ ] values over the study population at each HD stage. Multivariable estimation results using stage-specific, patient-specific and global approaches are presented.    Supplementary Material (Figures S3-S12). Tables 6 and 7 show a comparison between estimation errors obtained for the markers analyzed in this study and in previous studies. Tables S6-S9 in the Supplementary Material show mean absolute and root mean square errors for the analyzed markers.

Discussion
We investigated changes in QRS duration, amplitude and morphology at varying [K + ], [Ca 2+ ] and HR in ESRD patients. We designed [K + ] and [Ca 2+ ] estimators based on our proposed QRS morphological characteristics, taken both individually and in combination with T wave morphology markers. We showed the accuracy of our proposed estimators using three different approaches, stage-specific, patient-specific and global estimation, which outperformed previously proposed methods. Our results offered new non-invasive tools to monitor serum [K + ] and [Ca 2+ ], which could have a significant role in clinical practice and could contribute to reduce the mortality risk associated with abnormal electrolyte levels in ESRD patients.

Characterization of QRS Complex Amplitude, Duration and Morphology in ESRD Patients during and after HD
In ECG recordings of ESRD patients, we evaluated QRS duration and amplitude, measured by markers QRS w and QRS a , and QRS morphological characteristics, measured by markers d u w,Q , d a,Q , d NL w,Q and d NL a,Q , which was proposed here for the first time. All markers except for QRS w presented significant changes during and after HD. These changes were strongly associated with variations in [K + ] and [Ca 2+ ] but not in HR. The inconsistent relationship between QRS markers, such as QRS w , and HR has been investigated in previous works, including the study by Hnatkova et al. [46], who reported increases in QRS duration with increasing HR in 35% of their patients and decreases in QRS duration in the remaining 65%. In line with these results, we found that QRS became markedly wider at higher RR intervals for 34% of the patients; it became markedly narrower for 21% of the patients; its width moderately changed with RR for 21% of the patients; and it showed poor association with RR for the remaining 24% of the patients (see Table S2 in Supplementary Material). This led to a median r value of 0.42 between QRS w and RR over the 29 ESRD patients, reflecting a notably weaker relationship between QRS w and RR as compared to other depolarization markers. Regarding QRS amplitude, even if we found QRS a to remarkably change during HD, this marker depended on ECG amplitudes at specific time points, which, in noisy ambulatory recordings, could lead to large changes not associated with variations in electrolyte levels. On the other hand, changes in warping-based markers accounted for deviations in the whole QRS waveform and could thus be better suited for ambulatory monitoring. In particular, if QRS duration on top of amplitude changes occurred in the inter-dialytic interval, as, e.g., reported during advanced ischemia [47,48], these changes could be reflected in our proposed QRS warping markers.
On top of investigating the relationship between QRS markers and [K + ] or [Ca 2+ ], we also investigated the relationship between these markers and sodium concentration ([Na + ]).
[Na + ] variations were less remarkable than those of [K + ] and [Ca 2+ ] during HD in our patients (see Table S3 in Supplementary Material). Importantly, none of the markers were strongly associated with [Na + ] (see Figure S13 and Table S4 in Supplementary Material).
[Na + ] could have been expected to play a more important role in modulating depolarization markers because the fast sodium current was primarily responsible for phase 0 of the action potential and its changes might have manifested in QRS complex alterations. However, we found that the QRS complex became markedly wider at higher [Na + ] for 43% of the patients; it became markedly narrower for 7% of the patients; its width moderately changed with [Na + ] for 14% of the patients; and it showed poor association with [Na + ] for the remaining 36% of the patients (see Table S5 in Supplementary Material).
High inter-individual variability was found in the relationship between the analyzed QRS markers and [K + ] or [Ca 2+ ]. This was especially remarkable in the case of d a,Q , which presented high IQR in the intra-patient correlation coefficients with [K + ] and [Ca 2+ ]. Such high dispersion could be explained by QRS polarity effects, as a reduction (increase, respectively) in the absolute amplitude could be reflected as either positive or negative d a,Q , depending on QRS being predominantly positive or negative.
Changes in ECG characteristics induced by variations in [K + ] and [Ca 2+ ] have been extensively investigated in terms of ventricular repolarization. A number of studies have characterized changes in T wave width, amplitude, slope or slope-to-amplitude ratio [10,11,13,49,50]. We have recently quantified changes in T wave nonlinear dynamics and morphology and have shown their relationship with [K + ] and [Ca 2+ ] variations [28,29,31,33,51].
The analysis of electrolyte-induced alterations in ventricular depolarization remains, however, much more limited. In a study including 923 patients with severe hyperkalemia, sine wave-shaped QRS complexes were observed in almost 36.7% of patients [38]. In the present study, we could observe such behavior in most of the patients (Figure 4). Inconsistent results have been reported in relation to the effects of [K + ] on QRS duration, with a larger proportion of studies reporting QRS widening [1,34,[52][53][54][55][56] and others reporting QRS narrowing [2,35] with increased [K + ]. Here, we observed no significant changes in QRS width during and after HD. Regarding QRS amplitude, we found QRS a to be strongly negatively correlated with [K + ] and positively correlated with [Ca 2+ ], in accordance with the increase in QRS a with decreasing [K + ], as described by Astan et al. [35]. We observed similar results for the QRS morphology-based amplitude marker d a,Q . Our study characterized additional QRS morphological changes by the warping markers d u w,Q , d NL w,Q and d NL a,Q further extended these results to provide a robust characterization of QRS changes during and after HD in ESRD patients. Figure 8 illustrates QRS complexes and T waves at the start and end of HD.

Multivariable Predictors of [K + ] and [Ca 2+ ] Based on Depolarization and Repolarization Characteristics
Based on the novel QRS morphology markers of this study and on already proposed T wave markers [29,30,39] ] linear estimators, we used stage-specific, patient-specific and global estimation approaches. Overall, the stage-specific approach rendered results with both mean and median estimation errors very close to zero but with higher dispersion than in the patient-specific approach. The latter approach would be suitable for clinical application, as electrolyte levels could be predicted in each patient based on a short ECG recording of the same patient.
Multivariable estimators combining information from d u w,Q and d u w,T outperformed univariable estimators (the ones proposed by us and other authors [10][11][12]). The estimation errors for the multivariable estimators were lower (Figure 7, Figures S3-S12 in Supplementary Material, Tables 6 and 7) and the correlation coefficients r between actual and estimated electrolyte levels were higher than for univariable estimators, both for the patient-specific and global estimation approaches (see Tables 4 and 5 [29,30]. These results supported the use of our proposed QRS markers to improve prediction of [K + ] and [Ca 2+ ] by ECG repolarization markers. ECG depolarization-based estimators have been scarcely investigated in the literature for serum electrolyte monitoring. In [57], an ECG-based [K + ] estimator was designed using QRS duration in addition to T wave markers, but QRS duration was found not to be highly correlated with [K + ], in agreement with our present results for QRS w . Pilia et al. [14] reviewed studies evaluating QRS amplitude and width features, but no improved serum electrolyte prediction by incorporating these features into repolarization-based estimators was provided. Here, we proposed univariable and multivariable [K + ] and [Ca 2+ ] estimators that included information from the whole morphology of the QRS complex. By accounting for characteristics beyond QRS amplitude and width, these estimators could offer more robust performance for ambulatory monitoring of ESRD patients and overcome some limitations of previously proposed markers, such as their dependence on blood volume [36,37]. For all the estimators we built, we found that the values of [K + ] and [Ca 2+ ] at the beginning of each of the two HD sessions, i.e., time points h 0 and h 48 , were the more challenging to estimate, as could be observed from the largest estimation errors at those time points (Figure 4). Here, we duplicated the values of the ECG markers at h 0 and h 48 to give them more weight in the training step, as all other measures corresponding to h 1 , h 2 , h 3 and h − 4 were more similar to one another and the learning could otherwise be biased towards such measures.

Study Limitations and Future Research
This study investigated 29 ECG recordings of ESRD patients. Although the dataset was originally planned to include a larger number of patients, ECG acquisition had to be stopped due to the situation generated by the COVID-19 pandemic. Future studies should investigate the application of the proposed methods to larger numbers of patients.
For each patient, six blood samples were available, five of them taken during an HD session and the sixth one at the beginning of the following HD session. Future studies could be designed to have more frequent [K + ] and [Ca 2+ ] measurements, especially during the first HD hour, when electrolyte levels vary most remarkably. This could help to improve the learning of the estimators.
We used linear estimators due to the small number of samples, particularly when using a patient-specific approach. Future work could investigate the use of nonlinear estimators [11,32], which could prove to be particularly relevant for the estimation of [K + ] and [Ca 2+ ] at the start of HD sessions when these take values far from those at other HD stages.
We did not have access to measurements of blood volume or of other variables that could be used to infer them. Studies on other datasets where such measurements were available could test the relationship between the markers d u w,Q and d u w,T used in our [K + ] and [Ca 2+ ] estimators and the blood volume. In addition, some of the patients analyzed in the study had diseases, such as diabetes mellitus. We did not find significant differences in the analyzed markers between diabetic and non-diabetic patients. Nevertheless, future studies addressing larger patient cohorts could investigate the impact of diseases additional to ESRD on the relationship between QRS markers and electrolytes.
Some publications have reported the use of deep learning methods for hypo-and hyperkalemia screening from the ECG. While such methods require large amounts of training data and could not be addressed with the database analyzed here, further investigations could extend the present work to include machine learning techniques for the same purpose [58][59][60].
We focused our research on [K + ] and [Ca 2+ ] estimation. [Na + ] was found to present less notable variations during HD and none of our analyzed QRS markers showed significant association with it in the reduced dataset analyzed in this study, which should be further tested in larger patient cohorts. Although variations in other electrolytes, such as magnesium [Mg 2+ ], have also been shown to alter the ECG to some extent [2,24,[61][62][63], [Mg 2+ ] measurements were not available for the present study.
An additional future research line related to this work will be aimed at extending the present investigations to include patient-specific electrophysiological simulations in biventricular models embedded in torso models. From those simulations, realistic ECGs could be computed and used to gain understanding on the relationship between ECG characteristics and [K + ] or [Ca 2+ ].

Conclusions
Our proposed QRS morphology markers presented remarkable changes during and after HD, which were strongly associated with [K + ] and [Ca 2+ ] in the ESRD patients. Multivariable estimators based on combined QRS and T wave morphological variability allowed accurate prediction of [K + ] and [Ca 2+ ], outperforming estimators based on only ECG depolarization or repolarization. These results could pave the way to ambulatory, non-invasive monitoring of electrolyte levels, which could help to prevent fatal ventricular arrhythmias in ESRD patients.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/s22082951/s1, Figure S1: QRS complexes, including dominant and non-dominant polarity; Figure S2: Box plots of [K + ] and [Ca 2+ ] estimation errors e r ; Figures S3-S12: Bland-Altman plots between actual and estimated [K + ] and [Ca 2+ ]; Figure S13: Changes in sodium levels ([Na + ]) along HD stages; Table S1: Pearson correlation coefficient (r) to assess the association between QRS and T waves based multivariable estimators; Table S2: Pearson correlation coefficient (r) between QRS complex width and RR interval in 29 ESRD patients; Table S3: P-values to assess differences in [Na + ] between consecutive time stages; Table S4: Pearson (r) and Spearman (ρ) correlation coefficients between QRS complex markers and [Na + ] in 29 ESRD patients; Table S5: Pearson correlation coefficient (r) between QRS complex width and [Na + ] in 29 ESRD patients; Tables S6-S9: Mean absolute errors (MAE) and Root mean square errors (RMSE) using stage-specific (S), patient-specific (P) and global (G) approach based [K + ] and [Ca 2+ ] estimators. Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The dataset is still ongoing and it is available upon request to the corresponding author.