Nonlinear T-wave time warping-based sensing model for non-invasive personalised blood potassium monitoring in hemodialysis patients: a pilot study

sensors Article Nonlinear T-Wave Time Warping-Based Sensing Model for Non-Invasive Personalised Blood Potassium Monitoring in Hemodialysis Patients: A Pilot Study Flavio Palmieri 1,2,3 , Pedro Gomis 1,4 , José Esteban Ruiz 5, Dina Ferreira 3, Alba Martín-Yebra 2,6 , Esther Pueyo 2,6 , Juan Pablo Martínez 2,6 , Julia Ramírez 7and Pablo Laguna 2,6,*   Citation: Palmieri, F.; Gomis, P.; Ruiz, J.E.; Ferreira, D.; Martín-Yebra, A.; Pueyo, E.; Martínez, J.P.; Ramírez, J.; Laguna, P. Nonlinear T-Wave Time Warping-Based Sensing Model for Non-Invasive Personalised Blood Potassium Monitoring in Hemodialysis Patients: A Pilot Study. Sensors 2021,21, 2710. https://doi.org/10.3390/s21082710 Academic Editor: Piotr Augustyniak Received: 25 January 2021 Accepted: 30 March 2021 Published: 12 April 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). 1Centre de Recerca en Enginyeria Biomèdica, Universitat Politècnica de Catalunya, 08028 Barcelona, Spain; [email protected] (F.P.); pedr[email protected] (P.G.) 2CIBER en Bioingeniería, Biomateriales y Nanomedicina (CIBER-BBN), 50018 Zaragoza, Spain; [email protected] (A.M.-Y.); [email protected] (E.P.); [email protected] (J.P.M.) 3Laboratorios Rubió, Castellbisbal, 08755 Barcelona, Spain; dferr[email protected] 4Escuela Superior de Ingeniería, Ciencia y Tecnología, Universidad Internacional de Valencia, 46002 Valencia, Spain 5Nephrology Ward, Hospital Clínico Universitario Lozano Blesa, 50009 Zaragoza, Spain; [email protected] 6BSICoS Group, I3A, IIS Aragón, Universidad de Zaragoza, 50018 Zaragoza, Spain 7William Harvey Research Institute, Queen Mary University of London, London E1 4NS, UK; [email protected] *Correspondence: [email protected] Abstract: 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. Keywords: electrocardiogram; periodic component analysis; T-wave morphology; time warping; noninvasive potassium sensing; personalised medicine 1. 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 10to 20-fold higher than in ageand 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 Sensors 2021,21, 2710. https://doi.org/10.3390/s21082710 https://www.mdpi.com/journal/sensors Sensors 2021,21, 2710 2 of 13 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 – 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 – 17 ]. In previous studies [ 18 – 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 du w , dw , and ˆ dw,c (unsigned, signed, and heart rate corrected T-wave morphology variations in time, respectively), da (T-wave morphology variations in amplitude), and their non-linear components ( dNL w and dNL 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, πCB , or just restricting to the T-wave, πCT [ 20 ]. This work [ 20 ] showed that dw and ˆ dw,c had the highest correlation with ∆[K+] . They also showed that πCT presented higher robustness against noise than PCA or πCB , 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 2was 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 dw and ˆ dw,c (purple and green boxplots, respectively, in Figure 2of 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 πCT -based markers, dw and ˆ dw,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. Sensors 2021,21, 2710 3 of 13 Time (min) 5 0 ℎ0 2880 ℎ5 65 ℎ1 125 ℎ2 185 ℎ3 215/245 ℎ4 HD start HD end HD start min 5 ≡HD start on the morning of day 1 min 2880 ≡HD start on the morning of day 3 Figure 1. Data acquisition protocol. Holter electrocardiogram (ECG) signals of end-stage renal disease patients undergoing hemodialysis (ESRD-HD) patients were acquired throughout 48 h, starting 5 min before the beginning of the HD therapy. Six blood samples were collected at the beginning of the therapy ( h0 ), each hour during the HD ( h1 , h2 , h3 ), at the end ( h4 , at minute 215th or 245th, depending on the HD duration) and before the beginning of the next HD session (h5). Figure 2. Boxplots showing the distribution of ∆[K+] (blue) and the described πCT -based time warping biomarkers dw (purple) and ˆ dw,c (green), computed at each time points ( h0 to h5 ), see Figure 1. The central line of the boxplots represents the median, the edges of the box are the 25-th and 75-th percentiles, and the whiskers extend to the most extreme data points not considered as outliers. The notches represent the 95% confidence interval of the median, calculated as q2− 1.57 (q3−q1)/√n and q2 + 1.57 (q3−q1)/√n being q2 the median, q1 and q3 are the 25-th and 75-th percentiles, respectively, and nis the sample size. Finally, red “+” denotes outliers. Data adapted from [20,21]. 2. 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 Sensors 2021,21, 2710 4 of 13 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 . 3. 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. ECG filtering, QRS detection, T wave delineation 𝜋𝜋CA computation, QRS detection, T wave delineation Raw ECG (a) (b) (c) — — Time-warping analysis, 𝑑𝑑𝑤𝑤 and 󰆹 𝑑𝑑𝑤𝑤,𝑐𝑐computation [21] Windows selection and MWTW extraction Fitting models and � ∆𝑑𝑑,𝑚𝑚 𝑓𝑓𝐾𝐾+estimation (d) (e) (f) (g) Figure 3. Flow chart showing the ECG processing steps performed in this study. ( a ) Raw ECG (the eight independent leads I, II, V1 to V6 are shown) obtained from one of the enrolled ESRD-HD patients (see Section 2). ( b ) Preprocessed ECG as described in Section 3.1. ( c ) π CA is applied and both QRS complexes and T-waves (TW in the legend) are detected and delineated as detailed in Section 3.2. ( d ) From 2-min wide windows, ( e ) a mean warped T-wave (MWTW) is extracted and ( f ) T-wave morphology markers dw and ˆ dw,c are computed as stated in Section 3.3. ( g ) The fitting models for ˆ ∆f d,m[K+] estimation are evaluated as in Section 3.5. In this example, a cubic model with m=ais presented. Sensors 2021,21, 2710 5 of 13 3.1. 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]. 3.2. 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, πC 1, which was selected for subsequent analysis and T waves were again delineated by using the above-mentioned delineator [26]. 3.3. Warping-Based T-Wave Morphology Markers All T-waves from πC 1 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, dw and dw,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 ( h4 in Figure 1). A detailed description of the computation of the warping markers here analysed can be found in [ 21 , 23 ], describing how dw represents a relative measure of morphological changes between two T-waves. Likewise, ˆ dw,c is obtained from dw 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 ]. 3.4. 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 (h4), where a blood sample was taken: ∆[K+](hi) = [K+]hi−[K+]h4(1) being [K+]hi the concentration at the hi -th time point (see Figure 1) and [K+]h4 the reference concentration at the end of the HD treatment. The ∆[K+] distribution across patients for each hour is presented in Figure 2. 3.5. Marker Fitting Models for ∆[K+]Estimation For a given patient p , the relationship between the marker d∈ {dw , dw,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: ˆ ∆l d[K+](hi) = αld(hi), (2) ˆ ∆q d[K+](hi) = αqd(hi) + βqd2(hi), (3) ˆ ∆c d[K+](hi) = αcd(hi) + βcd2(hi) + γcd3(hi), (4) Sensors 2021,21, 2710 6 of 13 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 hi -th ∆[K+](hi) value from the training-set and evaluating the prediction error at this hi -th point, repeating this for all possible hiexclusions. 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 2in this paper and in Corsi et al. [16] in Figures 2 and 4 . That was implemented by imposing: ∂ˆ ∆d[K+] ∂d≥0, (5) which for positive values of the marker, d> 0, implies αl≥ 0, αq≥ 0, βq≥ 0, αc≥ 0, βc≥ 0, and γc≥ 0. 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. 3.6. Statistical Analysis Spearman’s rank and Pearson’s correlation coefficients ( ρ and r , respectively) were used for correlation analyses between ∆[K+] and ˆ ∆f d,m[K+] , where f∈ {l,q,c} denotes the fitting model, d∈ {dw,ˆ dw,c} the T-wave morphology parameter and m∈ {a,o} the estimation method. This analysis gives information about both the monotonic relation and the strength of the association between each modelled T-wave morphology parameter and ∆[K+] , thus providing a more complete characterisation. In addition, for each patient p and hour hi, an estimation error ef d,m(hi,p)was computed as: ef d,m(hi,p) =   ˆ ∆f d,m[K+](hi,p)−∆[K+](hi,p)  (6) where i∈ {0, 1, 2, 3, 5} is the set of hours where the computation of the estimation error is meaningful. Note that h4 is the reference point where both ∆[K+] and ˆ ∆f d,m[K+] are equal to zero and therefore it is excluded from error computation to avoid a biased error evaluation. The value ˆ ∆f d,a[K+](p , hi) represents the estimation at time hi when training is done including all available hi hours values from the corresponding patient, while ˆ ∆f d,o[K+](p , hi) represents the estimate at hour hi when all but the hi -th point of the patient are used in the training. All statistical analyses were performed using MATLAB version R2019a and results are given as the median and interquartile range (IQR). 4. Results The median and IQR values of intra-patient Spearman’s ( ρ ) and Pearson’s ( r ) correlation coefficients, computed between ∆[K+] , and ˆ ∆f d,m[K+] , are given in Table 1. The same table also displays the median and IQR values of the errors ef d,m(p , hi) , pooling together all patients and all blood extractions (ALL), and segregated for hour h0and h5. Sensors 2021,21, 2710 7 of 13 Table 1. Intra-patient ρ , r , ef d,m —either when pooling all patients and blood samples together (ALL) or specifically for h0 and h5 —evaluated between ∆[K+] and ˆ ∆f d,m[K+] , expressed as median (interquartile range (IQR)), for each model f∈ {l , q , c} , marker d∈ {dw,ˆ dw,c} , and estimation rule m∈ {a,o}. d f m ρref d,m ALL h0h5 dw la0.83 (0.33) 0.86 (0.35) 0.30 (0.48) 0.28 (0.77) 0.29 (0.55) o0.77 (0.48) 0.76 (0.47) 0.38 (0.61) 0.56 (1.10) 0.45 (0.66) qa0.83 (0.36) 0.91 (0.29) 0.22 (0.34) 0.24 (0.58) 0.27 (0.49) o0.83 (0.49) 0.77 (0.51) 0.38 (0.59) 0.64 (1.15) 0.63 (0.60) ca0.89 (0.35) 0.92 (0.27) 0.21 (0.34) 0.23 (0.37) 0.30 (0.54) o0.83 (0.49) 0.79 (0.61) 0.39 (0.72) 0.64 (1.24) 0.69 (0.75) ˆ dw,c la0.83 (0.31) 0.88 (0.34) 0.27 (0.50) 0.26 (1.03) 0.31 (0.54) o0.80 (0.44) 0.81 (0.34) 0.40 (0.63) 0.54 (1.11) 0.50 (0.59) qa0.83 (0.35) 0.90 (0.27) 0.21 (0.36) 0.25 (0.73) 0.27 (0.50) o0.80 (0.53) 0.77 (0.39) 0.41 (0.67) 0.57 (1.45) 0.71 (0.61) ca0.83 (0.31) 0.90 (0.25) 0.20 (0.39) 0.25 (0.67) 0.23 (0.52) o0.80 (0.49) 0.72 (0.45) 0.43 (0.81) 0.77 (1.25) 0.76 (0.80) Boxplots in Figure 4show the estimation error ef d,m(p , hi) distributions, sorted by hours hi , 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 h0 and h5 , 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 ( h4 ) 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 h0 ) and/or after (i.e., in case of h5 ) as opposed to h1 , h2 , and h3 ; 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 dw,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). Sensors 2021,21, 2710 8 of 13 Linear (a)d=dw Quadratic (b)d=dw Cubic (c)d=dw (d)d=ˆ dw,c(e)d=ˆ dw,c(f)d=ˆ dw,c Figure 4. Estimation error ( ef d,m(p , hi) ) distributions across patients for each hour hi and when pooling all samples together (ALL). Panels ( a , d ) show results for linear models f=l ; panels ( b , e ) show the quadratic and panels f=q ; and ( c , f ) show the cubic model f=c . Yellow dots represent individual error values when m = a , while light-blue ones denote those obtained when m = o . Corresponding boxplots are depicted on top of each distribution: The black ones represent the errors in m=awhile the red ones represents error in case of m=o. “+” denotes outliers. (a) (b) (c) Figure 5. Examples of cubic models (red dotted lines) computed for a given patient by imposing different parameter restrictions for leave-on-out cross-validation method, the corresponding equations are reported above each panel. The resulting model without restrictions on {αc , βc , γc} is in panel ( a ), while those from imposing αc≥ 0, or full constrained model are presented in ( b , c ) respectively. In each panel: The blue diamonds represent measured ∆[K+] values at the hours {h0 , h1 , h2 , h3 , h4 , h5} ; while red dots are the estimated ˆ ∆c dw,o[K+] corresponding to the computed dw used in the training set and computed at {h1 , h2 , h3 , h4 , h5} ; the green square is the estimated ˆ ∆c dw,o[K+] corresponding to the dw at h0 , the hour excluded from the training set in this example, and then the one with higher risk for error in the estimation. See that only full set of parameters forced to be positive result in a monotonic, physiologically plausible, function. 5. Discussion In this study, we analysed ECG signals from 29 ESRD-HD patients. We extracted Twave morphology parameters, dw and ˆ dw,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 dw and ˆ dw,c markers. The performance Sensors 2021,21, 2710 9 of 13 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=dw 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=dw 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=ˆ dw,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 dw and ˆ dw,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 dw with ˆ dw,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 ˆ dw,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, ˆ dw,c since it is already subjected to an heart rate correction estimation [21]. A reduction in the median and IQR estimation error for d=dw 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=ˆ dw,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 most distant hours from the reference point ( h0 and h5 in this work) are the most interesting and challenging for ∆[K+] estimation. Indeed, these two are the time points where the estimation errors are the highest and the error distributions are the widest, which is particularly true when they are taken out of the training set in the m = o case (red boxplot). That could be considered as an indication of the high uncertainty in predicting such values, especially when the values to be estimated do not have closer samples before and/or after, resulting in the wide IQR values reported for ef d,m at h0 and h5 , Table 1. In general, the IQR value for ef d,m at h0 decreases for dw/ˆ dw,c in m = a from 0.77/1.03 for linear to 0.58/0.73