Nonlinear T-Wave Time Warping-Based Sensing Model for Non-Invasive Personalised Blood Potassium Monitoring in Hemodialysis Patients: A Pilot Study

Background: End-stage renal disease patients undergoing hemodialysis (ESRD-HD) therapy are highly susceptible to malignant ventricular arrhythmias caused by undetected potassium concentration ([K+]) variations (Δ[K+]) out of normal ranges. Therefore, a reliable method for continuous, noninvasive monitoring of [K+] is crucial. The morphology of the T-wave in the electrocardiogram (ECG) reflects Δ[K+] and two time-warping-based T-wave morphological parameters, dw and its heart-rate corrected version dw,c, have been shown to reliably track Δ[K+] from the ECG. The aim of this study is to derive polynomial models relating dw and dw,c with Δ[K+], and to test their ability to reliably sense and quantify Δ[K+] values. Methods: 48-hour Holter ECGs and [K+] values from six blood samples were collected from 29 ESRD-HD patients. For every patient, dw and dw,c were computed, and linear, quadratic, and cubic fitting models were derived from them. Then, Spearman’s (ρ) and Pearson’s (r) correlation coefficients, and the estimation error (ed) between Δ[K+] and the corresponding model-estimated values (Δ^[K+]) were calculated. Results and Discussions: Nonlinear models were the most suitable for Δ[K+] estimation, rendering higher Pearson’s correlation (median 0.77 ≤r≤ 0.92) and smaller estimation error (median 0.20 ≤ed≤ 0.43) than the linear model (median 0.76 ≤r≤ 0.86 and 0.30 ≤ed≤ 0.40), even if similar Spearman’s ρ were found across models (median 0.77 ≤ρ≤ 0.83). Conclusion: Results support the use of nonlinear T-wave-based models as Δ[K+] sensors in ESRD-HD patients.


Introduction
Heart failure is among the most common cardiovascular complications in end-stage renal disease (ESRD) patients [1,2]. In hemodialysis (HD)-dependent ESRD (ESRD-HD) patients, the risk of cardiovascular mortality caused by electrical instability is 10-to 20-fold higher than in age-and gender-matched healthy subjects [3,4]. This remarkable association can be explained by the extreme fluctuations in blood potassium concentration ([K + ]) occurring in between HD sessions [5,6]. These changes in [K + ] are usually clinically silent and occur without warning to the patient or to the doctor in the absence of blood tests [7]. Therefore, continuous noninvasive monitoring of [K + ] variations (∆[K + ]) is of great importance [8] as it would provide risk warnings, improving an ever-growing clinical need.
The electrocardiogram (ECG) reflects the electrical activity of the heart in a noninvasive and inexpensive way. Electrocardiographic consequences of ∆[K + ] are well known [9][10][11]: The earliest effects appear as narrowed and peaked T waves [12], followed by changes in the QT interval duration (or its corrected version, QTc) [13] and in repolarisation complexity [14]. Several studies in the literature have attempted to estimate [K + ] through the analysis of T-wave morphology changes, quantified by features representative of the T-wave shapes, [15][16][17]. In previous studies [18][19][20][21][22], we proposed and investigated six T-wave morphological parameters quantifying T-wave morphology changes by means of time warping analysis [23] for continuous non-invasive ∆[K + ] monitoring. These six T-wave morphology parameters included d u w , d w , andd w,c (unsigned, signed, and heart rate corrected T-wave morphology variations in time, respectively), d a (T-wave morphology variations in amplitude), and their non-linear components (d NL w and d NL a ) as described in [21,23]. In addition, we tested two lead space reduction techniques [20], principal component analysis (PCA) and periodic component analysis (πCA) [24], this later implemented in two different versions: By exploiting the complete QRST complex periodicity, πC B , or just restricting to the T-wave, πC T [20]. This work [20] showed that d w andd w,c had the highest correlation with ∆[K + ]. They also showed that πC T presented higher robustness against noise than PCA or πC B , making it the most suitable lead space reduction technique for ∆[K + ] tracking during the HD session, as well as in the post therapy monitoring before the next HD session. Nevertheless, a quantitative relation between these T-wave morphological parameters derived from ECG analysis and ∆[K + ] has not yet been established for clinical use, which would allow a non invasive measurement of ∆[K + ] value.
The direct assessment of a marker as [K + ] surrogate by Pearson correlation analysis implies the assumption of a linear relation between them. However, previous works have reported that the reconstruction of [K + ] from the ECG significantly improves by employing a quadratic regression [16]. This result is compatible with the findings we reported in [20,21] where the same study population described in Section 2 was investigated according to the protocol in Figure 1. In Palmieri et al. [21], we observed a non-linear correspondence between ∆[K + ] and the T-wave time warping biomarkers d w andd w,c (purple and green boxplots, respectively, in Figure 2 of the present study). Therefore, we hypothesised that using patient-specific nonlinear models based on T-wave time warping-derived markers can provide better quantitative assessment of ∆[K + ]. The aim of this study is to derive and to evaluate nonlinear polynomial sensing models to estimate ∆[K + ] by using πC T -based markers, d w andd w,c . As a reference, a patient-specific linear model is also estimated for each marker.
This paper is organised as follows. First, we describe the study population, a database recorded during the interdialytic interval between two HD sessions. We next expand on the methodology to calculate the T-wave morphology parameters, as well as the proposed models to monitor [K + ] fluctuations from the ECG. Finally, we present and discuss the results, ending with conclusions and considerations for future research.  The notches represent the 95% confidence interval of the median, calculated as q 2 − 1.57(q 3 − q 1 )/ √ n and q 2 + 1.57(q 3 − q 1 )/ √ n being q 2 the median, q 1 and q 3 are the 25-th and 75-th percentiles, respectively, and n is the sample size. Finally, red "+" denotes outliers. Data adapted from [20,21].

Materials
The ESRD-HD study population included 29 patients from the Nephrology ward at Hospital Clínico Universitario Lozano Blesa (Zaragoza, Spain). Inclusion criteria were (i) 18-year-old or older patients, (ii) being diagnosed with ESRD, and (iii) undergoing HD at least three times per week, with venous or cannula access. The study protocol was approved by the Aragon's ethics committee (CEICA, ref. PI18/003) and all patients signed informed consent. All procedures and methods were performed in accordance with the Helsinki Declaration. Further details concerning the study protocol and clinical features of the study population can be found in [20,21]. A 48-h, standard 12-lead ECG Holter recording (H12+, Mortara Instruments, Milwaukee, WI, USA, sampling frequency of 1 kHz, amplitude resolution of 3.75 µV) was collected for each enrolled patient, with the acquisition starting 5 min before the HD onset and lasting until the next HD session, programmed 48 h later. Simultaneously, to determine [K + ], a total of 6 blood samples were collected just before starting the HD and every hour during the HD session (5 in total), with a last extraction immediately before the next HD session (Figure 1). The current number of patients included in the database is still limited, hence this work should be interpreted as an exploratory pilot study .

Methods
In this section, the different steps required for the processing of the ECG signals are described and summarised in the block diagram presented in Figure 3. Time-warping analysis, and ̂, computation [21] Windows selection and MWTW extraction Fitting models and � ∆ ,

ECG Pre-Processing
Baseline wander was removed with a 0.5-Hz cut-off high-pass filter, implemented with a forward-backward 6-th order Butterworth filter [25]. Residual noise out of the T-wave band was removed with a 40 Hz cut-off frequency forward-backward 3-th order low-pass Butterworth filter. QRS complexes were detected and T-waves delineated using a wavelet-based single-lead delineation method applied to each of the 12 leads [26].

Lead Transformation by Periodic Component Analysis, πCA
Periodic component analysis is a lead space reduction technique aiming to emphasise the periodic structure of a signal [24,27]. In this work, πCA was applied with a one-beat periodicity to maximise the T-wave beat-to-beat periodic components on the transformed signal, as explained in [20]. For each ECG recording, a transformation matrix Ψ πCA was estimated as detailed in [20], and applied to the 8 independent standard leads, obtaining a new set of 8 transformed leads, named periodic components. In this way and by ordering the transformed leads inversely to their associated eigenvalue, the most beat-to-beat periodic components appear projected onto the first component, πC1, which was selected for subsequent analysis and T waves were again delineated by using the above-mentioned delineator [26].

Warping-Based T-Wave Morphology Markers
All T-waves from πC1 were further low-pass filtered at 20 Hz using a forwardbackward 6-th Butterworth filter to remove remaining out-of-band frequency components. T-waves in 2-min wide windows centered around the 5-th minute and 35-th minute of each available hour were selected, from which a mean warped T-wave (MWTW) was computed from all T-waves in each window [21,23]. Finally, the two T-wave morphology parameters, d w and d w,c , were computed by comparing each MWTW with respect to a reference MWTW, selected at the end of the HD session, resulting in relative markers to the reference point at the end of HD (h 4 in Figure 1). A detailed description of the computation of the warping markers here analysed can be found in [21,23], describing how d w represents a relative measure of morphological changes between two T-waves. Likewise,d w,c is obtained from d w marker after being compensated for T-wave morphological changes not attributable to ∆[K + ] but to heart rate changes occurring between the reference and analysis points [20,21].

Blood Potassium Concentration Variations ∆[K + ]
The two proposed biomarkers, measured along time, have been associated with the corresponding relative variations in [K + ] with respect to the [K + ] at the reference point (h 4 ), where a blood sample was taken: being [K + ] h i the concentration at the h i -th time point (see Figure 1) and [K + ] h 4 the reference concentration at the end of the HD treatment. The ∆[K + ] distribution across patients for each hour is presented in Figure 2.

Marker Fitting Models for ∆[K + ] Estimation
For a given patient p, the relationship between the marker d ∈ {d w , d w,c } and ∆[K + ] measured along time was modelled by means of a linear (l), quadratic (q), and cubic (c) regression models for each patient to noninvasively calculate ∆[K + ] values, according to the following models: respectively. The coefficients α l , α q , β q , α c , β c , and γ c were estimated for each patient p and marker d by using a least square regression analysis between∆[K + ] and ∆[K + ] values. For each patient and marker, the parameters of the three models were estimated with two different approaches: (i) By using all the available ∆[K + ] values ("m = a") and (ii) by adopting a leave-one-out cross validation ("m = o") by excluding the h i -th ∆[K + ](h i ) value from the training-set and evaluating the prediction error at this h i -th point, repeating this for all possible h i exclusions. Finally, to avoid physiologically meaningless∆ d [K + ] trends, the three models in Equations (2)-(4) were computed with a constrained parameter estimation in order to guarantee a monotonically increasing relationship between∆[K + ] and d, as physiologically expected and corroborated by the marker trend evolution in Figure 2 in this paper and in Corsi et al. [16] in Figures 2 and 4 . That was implemented by imposing: which for positive values of the marker, The case with d < 0 is anecdotal, see Figure 2, and most likely is due to outliers, since they do not follow physiological interpretations of T-wave narrowing with increased potassium.

Statistical Analysis
Spearman's rank and Pearson's correlation coefficients (ρ and r, respectively) were used for correlation analyses between where i ∈ {0, 1, 2, 3, 5} is the set of hours where the computation of the estimation error is meaningful. Note that h 4 is the reference point where both ∆[K + ] and∆ All statistical analyses were performed using MATLAB version R2019a and results are given as the median and interquartile range (IQR).

Results
The median and IQR values of intra-patient Spearman's (ρ) and Pearson's (r) correlation coefficients, computed between ∆[K + ], and∆   Figure 4 show the estimation error e f d,m (p, h i ) distributions, sorted by hours h i , using the linear (Figure 4a,d), the quadratic (Figure 4b,e), and the cubic (Figure 4c,f) models. In addition, the aggregated distribution for all hours is presented with the label (ALL). The widest error distributions are obtained for hours h 0 and h 5 , whose median and IQR are given in Table 1. These time points are of great interest since: (i) The samples are the furthest from the reference (h 4 ) and (ii) when they are estimated by using the leave-one-out (m = o) approach they do not have any temporarily close samples before (i.e., in case of h 0 ) and/or after (i.e., in case of h 5 ) as opposed to h 1 , h 2 , and h 3 ; this together with the fact that their associated marker values are also the farthest from the rest, Figure 2. Therefore, it seemed worthy performing a detailed hour-based error analysis.
An example of cubic modeling results for a given patient with and without parameter constriction for monotonic∆ c d w ,o [K + ] behaviour with d is presented in Figure 5. Results are given for no restrictions on {α c , β c , γ c } ( Figure 5a); by just imposing α c ≥ 0 ( Figure 5b); and by the full constrained model (α c ≥ 0, β c ≥ 0, γ c ≥ 0) (Figure 5c).

Discussion
In this study, we analysed ECG signals from 29 ESRD-HD patients. We extracted Twave morphology parameters, d w andd w,c , previously reported to have a strong correlation with ∆[K + ] [20,21]. Then, we proposed and compared, the use of linear, quadratic, and cubic regression models for ∆[K + ] estimation from d w andd w,c markers. The performance of each model was evaluated through Spearman's and Pearson's correlation coefficients of the estimated∆[K + ] with respect to actual ∆[K + ] values and through hourly-based absolute estimation errors. The results on ESRD-HD patients here reported showed that non-linear regression models could be advantageously used to quantitatively estimate ∆[K + ] and could, therefore, be an effective tool for remote, frequent, and noninvasive monitoring of ESRD-HD patients.
According to Spearman's correlation coefficient (ρ) between measured and estimated variations in [K + ], similar ρ median values were found across the three models, with 0.06 being the highest median increment when moving from a linear to a cubic model for d = d w in m = a and, thus, denoting an analogous monotonic relationship between real ∆[K + ] and estimated values (∆ f d,m [K + ]). On the other hand, an improvement can be appreciated when comparing Pearson's correlation coefficient (r) evaluated in the three models, being the IQR reduced in d = d w and m = a by 0.06 and 0.08 when comparing the quadratic and cubic models, respectively, with respect to the linear model. Similar considerations can be made for d =d w,c . This is an expected outcome since the models here proposed were designed to avoid distorting the original monotonic increasing relationship between ∆[K + ] and the ECG derived markers, but only to adjust for the linear/non-linear relationship between them. However, the overall performance decreases considerably when the leave-one-out method, m = o, was used, being the median r lower and the IQR wider than in m = a. Also, for both d w andd w,c in m = o, a remarkable increase in the IQR can be observed when comparing linear and cubic models: From 0.47 to 0.61 for the first marker and from 0.34 to 0.45 for the second one. Overall, these findings seem to suggest that the cubic model does not provide any additional advantages to the linear or quadratic models in estimating ∆[K + ] using the leave-one-out approach. Therefore, the results we observed for m = a could potentially be affected by over-fitting.
Another interesting observation can be made when comparing d w withd w,c in terms of r: For the linear model and m = a, a small gain is obtained by heart rate correction, which is more significant for m = o. However, this improvement for the heart rate corrected index d w,c vanishes in m = a when using the quadratic model or the cubic model getting even worse in m = o. This can also be a result of the over-fitting in these estimates,d w,c since it is already subjected to an heart rate correction estimation [21].
A reduction in the median and IQR estimation error for d = d w in m = a results when hours and patients values are pooled together. The IQR decreases from 0.48 for the linear model to 0.34 for both the quadratic and cubic models. The median error goes from 0.30 in the linear model to 0.22 and 0.21 in the quadratic and the cubic models, respectively. An analogous trend can be found for d =d w,c in m = a: IQR reduces from 0.50 in f = l to 0.36 in f = q and to 0.39 in f = c. However, for both markers, the improvements disappear when the leave-one-out method m = o is used, which would support the previously hypothesised over-fitting for m = a. These outcomes would point at the quadratic model as the most suitable model for ∆[K + ] estimation in m = a, as well as in m = o, even if in this latter case the advantage is not very remarkable. Moreover, as mentioned above, there is no clear benefit in using a cubic rather than a quadratic model in any of both m = a and m = o cases, probably due to the full constrained parameter estimation rule we imposed which, when applied to the cubic model, we observed it resulted in a very small cubic term, reducing to quadratic model as in Figure 5. The results observed so far may lead to the conclusion that, according to the performance metrics r or e f d,o considered, the observed improvement for quadratic model estimation in the case of m = a vanishes, or it is largely attenuated, in m = o. However, when analysing data distributions we realise that values of d w andd w,c markers are not evenly distributed in all the analysed range (see Figure 2). This fact can imply an overweight of small d values in m = o modelling, penalising the estimates at h 0 and h 5 , which present d values that might not be well represented in the training set. This could also mean that the leave-one-out cross-validation needs to be cautiously framed when the value of d to be estimated is far from those used in the training set range, which in our dataset usually happens at h 0 and/or at h 5 as exemplified in Figure 6. In these cases, the estimation error between real ∆[K + ] and∆ In the following, some limitations of our study are acknowledged. If blood samples had been collected more frequently during the early stage of the HD treatment when [K + ] and, consequently, d more rapidly change-covering a broad range of values-then the model training set in m = o could have better represented all the possible cases of d in the quadratic as well as in the cubic model, and then the results could have been more conclusive for the non-linear modelling improvement in predicting [K + ]. If this refined learning would have been done, or is done in future studies, it will, predictably, result in less error at the extreme times h 0 and h 5 of the process, and consequently also in a notably improved performance of the quadratic model both for m = a and for m = o.
Another limitation that should be taken into account when interpreting this work's results is the lack of perfect time synchronisation between the actual ∆[K + ] and the evaluated d used for estimation at h 5 . As previously reported in [21], 44 h is the average ECG duration in our database-not 48 h, when the last blood sample is takenl-mainly due to electrode detachment or early battery exhaustion. However, in a recent study [20], we observed a low marker dynamics in the late post-HD treatment. Therefore, with some degree of confidence, we have assumed that the estimation error obtained between ∆[K + ] and∆ f d,m [K + ] at h 5 would be quite similar if the actual value-had the ECG lasted, as planned, for 48 h-had been used for modelling.
Specific aspects of ESRD-HD patients' clinical status (e.g., possibility of previous infarction not always revealed in clinical history) could have influenced the results, generating the inter-patient variability here observed. In addition, the accuracy of the proposed models in estimating potassium variations for patients other than ESRD-HD remains to be assessed.
Finally, the reduced number of patients and available blood samples for each patient included in this study also represents a limitation to frame the conclusion of the work. Indeed, even if the proposed approach may entail a significant step towards a robust and reliable ∆[K + ] sensing from time-warping based biomarkers, it needs to be validated in larger cohorts before any translation to clinical practice. However, the available data would suggest that a patient-specific quadratic model could estimate ∆[K + ] time trends with better accuracy than a linear-model. Also, in real practice, this method implies the collection of several blood samples, which may result in cumbersome procedures. It remains to be studied to what extent the models learned in one session can be extrapolated for sessions in later days/weeks, reducing the learning to just a single session.
Future studies should be conducted in a larger population including not only ESRD-HD patients but also subjects at risk of [K + ] imbalance, such as those with diabetes mellitus [28] or severe cardiovascular events like myocardial infarction [29]. In addition, the proposed estimation models should be validated in a follow-up study where the models are learned at the initial HD session and used in later HD sessions to measure ∆[K + ]. In such studies the complete learning with m = a at the initial HD session could be evaluated by its prediction value at subsequent sessions, without any overfitting risk. At this future analysis, we expect that m = a approach will show better performance, in terms of correlation and estimation error, than the one reported here for the m = o case, since the models' coefficients will be estimated over the six ∆[K + ] values (and not just over five as in m = o), thus covering the full range of d values for each patient.

Conclusions
The present study showed the advantage in using non-linear models in estimating ∆[K + ] measurements in ESRD-HD patients based on T-wave-derived markers. These results suggest a new noninvasive strategy for ECG-based [K + ] sensing, with large implications for monitoring patients with cardiovascular and renal diseases, providing a meaningful tool for a personalised ambulatory cardiac risk assessment of these patients.

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