Monitoring blood potassium concentration in hemodialysis patients by quantifying T-wave morphology dynamics
Abstract
We investigated the ability of time-warping-based ECG-derived markers of T-wave morphology changes in time (dw) and amplitude (da), as well as their non-linear components (dNLw and dNLa), and the heart rate corrected counterpart (dw,c), to monitor potassium concentration ([K+]) changes (Δ[K+]) in end-stage renal disease (ESRD) patients undergoing hemodialysis (HD). We compared the performance of the proposed time-warping markers, together with other previously proposed [K+] markers, such as T-wave width (Tw) and T-wave slope-to-amplitude ratio (TS/A), when computed from standard ECG leads as well as from principal component analysis (PCA)-based leads. 48-hour ECG recordings and a set of hourly-collected blood samples from 29 ESRD-HD patients were acquired. Values of dw, da, dNLw, dNLa and dw,c were calculated by comparing the morphology of the mean warped T-waves (MWTWs) derived at each hour along the HD with that from a reference MWTW, measured at the end of the HD. From the same MWTWs Tw and TS/A were also extracted. Similarly, Δ[K+] was calculated as the difference between the [K+] values at each hour and the [K+] reference level at the end of the HD session. We found that dw and dw,c showed higher correlation coefficients with Δ[K+] than TS/A—Spearman’s (ρ) and Pearson’s (r)—and Tw—Spearman’s (ρ)—in both SL and PCA approaches being the intra-patient median ρ≥0.82 and r≥0.87 in SL and ρ≥0.82 and r≥0.89 in PCA respectively. Our findings would point at dw and dw,c as the most suitable surrogate of Δ[K+], suggesting that they could be potentially useful for non-invasive monitoring of ESRD-HD patients in hospital, as well as in ambulatory settings. Therefore, the tracking of T-wave morphology variations by means of time-warping analysis could improve continuous and remote [K+] monitoring of ESRD-HD patients and flagging risk of [K+]-related cardiovascular events. Palmieri, F.; Gomis, P.; Ferreira, D.; Ruiz, J.E.; Bergasa, B.; Martín-Yebra, A.; Bukhari, H.A.; Pueyo, E.; Martínez, J.P.; Ramírez, J.; Laguna, P.
Full text
ͷ ǤǣȋͰͱͲͳʹ͵Ȍ Ƥ | (2021) 11:3883 | ǣȀȀǤȀͷͶǤͷͶ;ȀͺͷͻͿ;ǦͶͷǦ;ͿͻǦͻ www.nature.com/scientificreports Ǧ ͷǡǡ*ǡͷǡͺǡǡ±ͻǡͻǡ ÀǦǡͼǡǤǡͼǡǡͼǡÀǡͼǡ ÀͽƬǡͼ ǦǦǦǦ ȋ d w Ȍȋ daȌǡǦȋ dNL wdNL aȌ ǡ ȋ dw,cȌǡȋ [ K+ ] Ȍ ȋ [ K+ ] ȌǦȋȌȋȌǤ Ǧǡ [ K+ ] ǡǦȋ Tw ȌǦǦǦȋ T S/ A Ȍǡ ȋȌǦǤͺ;Ǧ ǦͿǦǤ dw ǡ daǡ dNL wǡ dNL adw,c ǦȋȌǡ Ǥ Tw T S/ A Ǥǡ [ K+ ] ơ [ K+ ] [ K+ ] Ǥ dw dw,cƥ [ K+ ] T S/ A Ȅǯȋ ρ ȌǯȋrȌȄ Tw Ȅǯȋ ρ ȌȄ Ǧ ρ ≥ 0 . 82 r ≥ 0 . 8 7 ρ ≥ 0 . 82 r ≥ 0 . 89 ǤƤ dw dw,c [ K+ ] ǡ ǦǦ ǡǤǡǦ Ǧ [ K+ ] Ǧ ƪ [ K+ ] ǦǤ Chronic kidney disease (CKD) is defined as the presence of kidney damage, persisting for 3 months or more, irrespective of the cause1. It represents a state of progressive loss of kidney function ultimately resulting in need for renal replacement therapy such as hemodialysis (HD) or transplantation. The development of CKD and its progression to this terminal stage, called end-stage renal disease (ESRD), remains a significant source of reduced quality of life and premature mortality2. In particular, sudden cardiac death (SCD) represents an important cause of death in ESRD-HD patients3. Various risk factors may be responsible for SCD in this patient population, including left ventricular hypertrophy and fibrosis, disordered bone-mineral metabolism, HD-induced changes in electrolyte, and fluid and acid-base status, which may lead to electrocardiographic (ECG) abnormalities and ventricular arrhythmia3,4. Recent studies have shown that blood potassium concentrations ( [ K+ ] ) outside the physiological interval are associated with increased mortality risk5. In healthy conditions, the maintenance of [ K+ ] homeostasis is ensured by normal renal activity6. However, ESRD-HD patients suffer from [ K+ ] imbalance, leading to a high incidence ͷCentre de Recerca en Enginyeria Biomèdica, Universitat Politècnica de Catalunya, Barcelona, Spain. CIBER en ÀǡȋǦȌǡǡǤLaboratorios Rubió, Castellbisbal, Barcelona, Spain. ͺValencian International University, Valencia, Spain. ͻNephrology Department, Hospital ÀǡǡǤͼǡ ǡ×ǡǡ ǡ Ǥ ͽWilliam Harvey Research Institute, Queen Mary University of London, London, UK. *email: ƪǤǤ
Vol:.(1234567890) Ƥ | (2021) 11:3883 | ǣȀȀǤȀͷͶǤͷͶ;ȀͺͷͻͿ;ǦͶͷǦ;ͿͻǦͻ www.nature.com/scientificreports/ of arrhythmic events. The pro-arrhythmic consequences of [ K+ ] imbalance can be explained considering that potassium currents are involved in the repolarization process of the cardiac action potential (AP), determining membrane potential and refractoriness of the myocardium7. Therefore, even modest deviations of [ K+ ] from its normal range (hypokalemia if [ K+ ] < 3.5 mmol/L or hyperkalemia if [ K+ ] > 5 mmol/L) may lead to hospitalisation or death in ESRD-HD patients8. Evaluation of [ K+ ] levels is currently based on blood samples that require further analyses in the laboratory, limiting continuous monitoring. Non-invasive markers able to track variations in [ K+ ] levels are therefore needed. The electrocardiogram (ECG) is a non-invasive, easily accessible, and inexpensive practice that reflects the electrical activity of the heart. In particular, the T-wave reflects the spatio-temporal repolarization of the ventricle, and its analysis has been used to measure the vulnerability of a patient to ventricular arrhythmias9. This fact is of particular interest because T waves are frequently altered in ESRD-HD patients4. The QT interval is the standard index of ventricular repolarization, and it has been proposed to monitor ESRD-HD patients10. However, the effects of HD on QT interval, and its corrected version QTc, are still controversial, since several studies11 reported a prolongation during the HD sessions, but others reported opposite trend or even no changes at all12. This motivates the analysis of the overall T-wave morphology as a potential potassium level marker. Different T-wave morphology markers have been previously reported to be correlated with [ K+ ] , such as the T-wave right slope13, the width of the T-wave ( Tw )14, the T-wave slope-to-amplitude ratio ( T S / A)15, and a morphology combination score, which integrates features like T-wave asymmetry, flatness and notching16. However, these markers rely on specific local features of the T-wave rather than in the overall T-wave morphology, which may have a stronger potential in following [K+ ] than indices based on local features. A recent study reported a time-warping based methodology to quantify changes in the overall T-wave morphology17. Six indices were proposed, du w and d a , reflecting morphological variations in time and amplitude, respectively, as well as their non-linear version, dNL w and dNL a as reported in17 and two novel markers derived from du w and named d w and dw , c . The main goal of this study is to investigate the potential of these markers in monitoring both hypoand hyperkalemia events excluding the variability due to the heart rate (HR) and to compare their performance against Tw and T S / A in standard single-lead approach and by applying principal component analysis (PCA) as multilead space reduction technique. However, some of the above mentioned indices may not be robust enough for our purpose. It is the case of d w which does not provide information about the direction of the T-wave morphological variation (i.e. if there is stretching or shortening) and has been found to be correlated with HR. Therefore, we have adapted the original methodology17 to account for hypoand hyperkalemia, and we propose a new marker that is independent of HR, thus offering a more precise [ K+ ] monitoring tool for arrhythmic risk stratification in ESRD-HD patients. Preliminary results extracted from a smaller subset of patients have been presented at Computing in Cardiology conference18,19 while the electrophysiological basis was studied in Bukhari etal.20. The novelties of the present study with respect to the state-of-the-art are: (1) the usage of T-wave time warping analysis for non-invasive [ K+ ] monitoring, together with the development of a HR correction tool for the time-warping marker, dw , c ; (2) the proposal of a PCA spatial transformation lead for marker extraction and (3) the validation of the proposed markers in comparison with previously published biomarkers ( Tw and T S / A ) and with their extraction from standard leads. Ǥ The study population included 29 patients from the Nephrology ward from Hospital Clínico Universitario Lozano Blesa (Zaragoza, Spain). Inclusion criteria were (i) 18-year-old (or older), (ii) having a diagnosed ESRD pathology and (iii) undergoing HD at least three times per week (with venous or cannula access). Table1 shows the population characteristics. The study protocol was approved by the Aragon’s research ethics committee (CEICA, ref. PI18/003) and all patients and/or their legal guardians signed informed consent. All the procedures and all the methods were performed in accordance with the Helsinki Declaration. The database collection is still ongoing, with the current size significant enough for a pilot study21,22. Ǥ General information. Sex, age, concomitant therapies (e.g. assumption of anti-arrhythmic drugs), kidney disease etiology and HD treatment related information were collected for each enrolled patient, as detailed in Table1. Blood sample analysis. For each patient, six blood samples were taken and analysed during the HD session: the first one at the HD onset and the next three, every subsequent hour (Fig.1, h 0 to h 3 in red). The 5-th blood sample was collected at the end of the HD (minute 215-th or 245-th, depending on the HD session duration) while the 6-th blood sample was taken after 48 h, immediately before the next HD session. Potassium, magnesium, calcium, urea, creatinine, bicarbonate and pH were measured from each blood test. Blood potassium concentrations values for each blood test are given in Table2. ECG measurements. A 48h, standard 12-lead ECG Holter recording, (H12+, Mortara Instruments, Milwaukee, WI, USA, sampling frequency of 1 kHz, amplitude resolution of 3.75 μ V), was obtained for each enrolled patient, starting the acquisition 5 min before the HD onset (Fig.1, blue line). The block diagram presented in Fig.2a describes the main steps of the whole data processing implemented in this work. ǦǤ ECG filtering. Holter ECG signals contain baseline drift and other noises, such as power-line and muscular activity (Fig.2a). Therefore, an initial pre-processing is needed to improve the signal-
Vol.:(0123456789) Ƥ | (2021) 11:3883 | ǣȀȀǤȀͷͶǤͷͶ;ȀͺͷͻͿ;ǦͶͷǦ;ͿͻǦͻ www.nature.com/scientificreports/ to-noise ratio (SNR) and enable ECG waveform analysis. First, baseline wander was removed with a high-pass, forward-backward 6-th order Butterworth filter with 0.5 Hz cut-off frequency23 (Fig.2a). Then, residual noise out of the T-wave band was removed with a 6-th order low-pass Butterworth filter with 40 Hz cut-off frequency. ECG waveform detection and delineation. A wavelet-based single-lead method24 was applied to detect QRS complexes and then delineate T-wave onsets and ends in each of the 12 leads. The wavelet transform (WT) decomposes the signal in the time-scale domain, allowing its representation at different resolutions. It is, therefore, a suitable tool to analyze ECG signals, which contain patterns with different frequency content (QRS complexes, P and T-waves). Single-lead delineation. The discrete dyadic WT is implemented in such a way that it keeps temporal resolution at different scales. The detection of the fiducial points is carried out across the adequate WT scales, attending to the dominant frequency components of each ECG wave: Q,R,S waves correspond to a simultaneous effect in scales 21–22 , while the T and P waves affect mainly scales 2 4 or 25 , see24 for details. ECG wave peaks correspond to zero crossings in the WT, and ECG maximum slopes correspond to WT’s maxima and minima. Depending on the number and polarity of the slopes found, a wave morphology is assigned and boundaries are located using threshold-based criteria. The onset (end) of a wave occurs before (after) the first (last) significant slope associated with the wave24. Selection rules for multi-lead delineation. To obtain multilead peak locations, a median post-processing selection rule over the single-lead-based detected locations is used. The post-processing rules for boundaries consist of ordering the single-lead annotations and selecting as the onset (end) of a wave the first (last) annotation whose k nearest neighbours lay within a δ ms interval24,25. Single-lead analysis. First, we performed the analysis using the single-lead ECG, taking the T-waves from leads V3 to V6, as used in a previous study26 for [ K+ ] estimation, and lead II being the most widely used in patient monitoring27. These T waves were further delineated by using the above motioned delineator24 and the biomarker estimation is performed as described below in section named “Time warping analysis”. Spatial lead reduction by principal component analysis. Next, a spatial lead reduction by Principal Component Analysis (PCA) was made since it was found to be a robust spatial transformation to emphasize waveform SNR28. In this work PCA was spatially applied to the 8 independent leads, learned over the T-wave segment to mainly emphasize this waveform, and resulting in 8 principal components (PCs) or transformed leads. The coefficients defining the PCA transformation were obtained from the eigenvectors of the 8× 8 interlead auto-correlation matrix computed over the T-waves in a 10-min wide window at the end of the HD session. The correct delineation of T-waves is crucial to emphasize only T-wave energy content. The first PCA, denoted as PC1, was used for the subsequent ECG analysis, as it is the transformed lead where the T-waves have maximal energy, and thus, maximal SNR for morphological characterisation28,29. PC1 was further delineated by applying again24, and each T-wave was further low-pass filtered at 20 Hz using a 12-th order Butterworth filter to restrict shape analysis to the dominant band of the T-wave so removing remaining noisy components that could still corrupt the T-wave shape analysis. ǦǤ Two-minute ECG segments, centred on the 5-th and 35-th minutes of each available hour, were analysed. The window duration was short enough to hold the assumption of stability for both [ K+ ] and HR values. Figure5a shows the average RR interval for each selected i-th 2-min segments for a given patient along the ECG recording. While the blood samples (purple diamonds) were collected each hour during the HD, the warping parameters were computed every half an hour to get a more detailed view over time. For each i-th 2-min segment, a mean warped T-wave (MWTW) was computed. First of all, the predominant T-wave polarity (e.g. upward, downward etc) within a given window, was defined as that having the highest number of occurrences. This polarity change can be physiological or induced by delineator oscillation when byphasic to regular T-waves appears almost indistinguishable. A T-wave was considered to have inverted polarity if the magnitude of its peak had negative sign and vice-versa. Only those T-waves having the same polarity as ECG acquisition 0565 125 185 215 245 HD Post HD 2880 ℎ1 ℎ0ℎ2ℎ3ℎ4ℎ5 Time (min) Figure1. Diagram of the study protocol: h 0 to h 5 are the time points (in minutes) for blood sample extraction. h 4 is taken at the end of the HD (minute 215-th or 245-th, depending on the HD duration).
ͺ Vol:.(1234567890) Ƥ | (2021) 11:3883 | ǣȀȀǤȀͷͶǤͷͶ;ȀͺͷͻͿ;ǦͶͷǦ;ͿͻǦͻ www.nature.com/scientificreports/ Original ECG ECG filtration PCA Time-warping analysisECG pre-processing (a) (b) Figure2. Analysis stages performed in this study. In panel(a) is the flow chart showing the ECG processing steps for T-wave time-warping markers extraction. The analysis starts with the original ECG, followed by a filtering step before spatial PCA analysis, to conclude with markers computation. Panel(b) shows an example of the linear and nonlinear time-warping markers for the same patient as in Fig.5a. In particular, subpanel (i) shows both the reference (blue) and the i-th MWTW (red) while subpanel (ii) shows the warping function (red dotted line) that optimally relates the reference and studied MWTWs. Subpanel (iii) shows the MWTWs after warping and subpanel (iv) are the normalized reference and warped MWTWs.
ͻ Vol.:(0123456789) Ƥ | (2021) 11:3883 | ǣȀȀǤȀͷͶǤͷͶ;ȀͺͷͻͿ;ǦͶͷǦ;ͿͻǦͻ www.nature.com/scientificreports/ the predominant one were considered in the following steps. Then, all these selected T-waves were aligned with respect to their gravity center and used to compute an initial MWTW 17. Finally, all the T-waves were checked to find and discard the outliers, defined as having T-wave duration outside the range T dmi±1.5 ×σdi centered at the i-th ensemble T-wave mean, T dmi , and bounded function of the T-wave duration standard deviation σdi . Among the remaining, only those T-waves highly correlated (Pearson’s correlation coefficient > 0.98) with the previous initial MWTW were used to recalculate the final MWTW. The MWTW at the end of the HD treatment was taken as the reference, given that it is the time when the patient (a) is supposed to have recovered the normal [ K+ ] level and (b) was discharged from hospital, being an appropriate reference for out-of-hospital ambulatory monitoring. Since hyperkalemia has been reported to cause T-wave inversions30, any MWTWs with negative-polarity was inverted before performing the warping with the reference MWTW. Previous to warping, the two MWTWs were aligned with respect to their gravity center, so that only changes in the T-wave morphology, and not those associated with their relative delay, were quantified by the warping algorithm. For comparison purposes, both Tw 14 and T S / A15 were extracted from each MWTW and their performance, with respect to T-wave time-warping based biomarkers in monitoring [ K+ ] , was assessed. This work perform a clinical study following previous analysis testing the marker by electrophysiological simulations as reported in Bukhari etal.20. Table 1. Characteristics of the study population. Values are expressed as number ( % ) for categorical variables, and median (IQR) for continuous variables. (N = 29) Age (years) 7 5 ( 12 ) Gender (male) 20 ( 70% ) Anti-arrhythmic drugs (yes) 9 ( 31% ) Implanted pace-maker (yes) 1 ( 3% ) Time under HD treatment (months) 15 ( 59 ) HD session duration 210 min 3 ( 10% ) 240 min 26 ( 90% ) Kidney disease etiology Diabetes mellitus 17 ( 59% ) Interstitial nephritis 2 ( 7% ) Glomerulonephritis 2 ( 7% ) Tuberous sclerosis 1 ( 3% ) Polycystic kidney 1 ( 3% ) Cancer 1 ( 3% ) Unknown 5 ( 18% ) HD liquid composition Potassium (1.5 mmol/L) 21 ( 72% ) Potassium (3 mmol/L) 5 ( 17% ) Potassium (decreasing) 3 ( 11% ) Calcium (2.5 mg/dL) 21 ( 72% ) Calcium (3 mg/dL) 8 ( 28% ) HD techniques Conventional 18 ( 62% ) Online 8 ( 28% ) Acetate-free biofiltration with decreasing intra-HD [K+]3 ( 10% ) Table 2. Blood potassium concentration [ K+ ] values (in mmol/L) at each blood extraction during the HD ( h 0 to h 4 ) and HR (beats/min). Spearman’s ( ρ ) and Pearson’s (r) intra-patient correlation coefficients between [ K+ ] and RR. Values are expressed as median (IQR). h 0 h 1 h 2 h3h 4 ρ r (K + )5.0 (1.4) 3.8 (1.1) 3.6 (0.8) 3.4 (0.7) 3.3 (0.6) 0.10 (1.35) 0.09 (1.45) HR 81 (28) 76 (28) 80 (23) 80 (17) 80 (25)
ͼ Vol:.(1234567890) Ƥ | (2021) 11:3883 | ǣȀȀǤȀͷͶǤͷͶ;ȀͺͷͻͿ;ǦͶͷǦ;ͿͻǦͻ www.nature.com/scientificreports/ T-wave time warping. The method here applied was originally proposed by Ramírez et al.17. Let fi(ti)=[fi(ti(1)),...,fi(ti(Ni))] T be the MWTW of a given i-th segment, and fr(tr)=[fr(tr(1),...,fr(tr(Nr))] T the reference MWTW, where t i=[ti(1),...,ti(Ni)] T and tr=[tr(1),...,tr(Nr)] T with N i and Nr being the total T-wave duration, in samples, of ti and tr respectively. Figure2b illustrates the warping method applied between one of the i-th MWTW (red) and the reference MWTW (blue). Let γ i(tr ) be the warping function that relates tr and ti , such that the composition ( fi◦ γ i)(tr) denotes the re-parametrization or time domain warping of fi(ti) using γ i(tr ) , i.e. ( fi◦ γ i)(tr) represents the amplitude values of fi(ti) if its temporal vector was tr . The square-root slope function (SRSF) was proposed instead of the original T-waves31,32 to find the optimal warping function. This was applied by performing time-warping on the SRSFs of the T-waves, preventing the “pinching effect” in cases when T-wave amplitudes differ33. This transformation is defined as: The optimal warping function is the one that minimizes the amplitude difference between the SRSF of fr(tr ) and fi( γ i(tr) ) 32: The dynamic programming algorithm was used to obtain the solution of this optimisation problem34. Figure2b(ii) shows the optimal warping function between the two waves in Fig.2b(i). The warped T-wave, fi(γ ∗ i (tr) ) is shown in Fig.2b(iii), together with the reference T-wave, fr(tr ) . Time warping biomarkers. The index du w (corresponding to the index denoted as d w in17), shown as the yellow area in Fig.2b(ii), quantifies the amount of warping needed to optimally align the two T-waves, and is defined as the average of the absolute difference value between γ∗ i (tr ) and tr: The original definition of du w(i ) 17 was modified here to allow the marker to be signed, therefore distinguishing T-wave widenings from narrowings. This signed dw(i ) was defined as: where sd ( i ) was used to account for the sign of the dw(i ) and it was computed as: with Nu r being the set of T-wave up-slope samples. A positive sign means that the fi(ti) has to be widened to fit the fr(tr ) and vice-versa for a negative sign. After applying time warping between both MWTWs, the amplitude difference between fr(tr ) and fi(γ ∗ i (tr) ) is quantified as the area contained between fr(tr ) and fi(γ ∗ i (tr) ) , normalized by the L2-norm of fr(tr ) : where sa(i)= N r n =1(fi(γ ∗ i (tr)) −fr(tr) ) is used to account for the da(i ) sign estimation. Both dw(i ) and da(i ) incorporate information from the linear and non-linear differences between both T-waves in time and amplitude domain, respectively. The non-linear components can be quantified as in17: where γ ∗ i , l(tr ) (green line in Fig.2b(ii)) is the best linear fitting to γ ∗ i (tr ) according to the least absolute residual criterion35. The parameter dNL w(i ) quantifies the non-linear warping by computing the area of the dashed magenta region between γ ∗(tr ) and γ ∗ i , l(tr ) (in Fig.2b(ii)). Finally, the marker dNL a(i ) quantifies the residual information in amplitude domain after normalising MWTWs (Fig.2b(iv)). (1) qf(t)=sign˙ f(t) ˙ f(t) . (2) γ∗ itr=arg min γi(tr) qfrtr−q[fi◦γi]tr =arg min γ i(tr) qfrtr−qfiγitr˙γi(tr) . (3) d u w(i)= 1 Nr Nr n=1 |γ∗ i(tr(n)) −tr(n)|. (4) dw(i)=sd(i) |sd(i)|1 Nr Nr n=1 |γ∗ itr(n)−tr(n)| . (5) sd(i)= n∈Nu r (γ ∗ itr(n)−tr(n))+ n/∈Nu r (tr(n)−γ∗ itr(n)). (6) da(i)= sa(i) sa(i) fi(γ ∗ i(tr)) −fr(tr) fr(tr) ×100 . (7) dNL w(i)=1 Nr Nr n=1 |γ∗ i(tr(n)) −γ∗ i,l(tr(n))| . (8) dNL a(i)= fr(tr) fr(tr)− fi(γ ∗ i(tr)) fi(γ ∗ i (tr)) ×100 .
ͽ Vol.:(0123456789) Ƥ | (2021) 11:3883 | ǣȀȀǤȀͷͶǤͷͶ;ȀͺͷͻͿ;ǦͶͷǦ;ͿͻǦͻ www.nature.com/scientificreports/ Heart-rate-corrected T-wave warping. It is well known that T-wave duration and QT interval are strongly dependent on HR36. Although aligning the T-waves according to their gravity centre reduces most of the dependence of dw(i) on HR, there may still be some residual dependence in T-wave morphology that should be compensated for (e.g. see Fig.5a around hours h = 9, 12 and 43). We assume that dw(i ) , as originally proposed in (4), can be modelled as the sum of two components: where dw , HR(i ) is the HR dependent component and dw , c(i ) is the non-HR dependent component accounting for ( K+ ) induced variations and possibly others not HR related. To estimate the corrected component dw , c(i ) we depart from the literature, where several formulae for HRdependency correction of repolarization related time intervals, like the QT interval, have been developed37–40, including a variety of approaches (e.g. linear, hyperbolic, exponential models etc.) being investigated and tested in view of the complex relationship between QT interval and HR39. To derive a correction formula and estimate dw , c(i ) , we started from a linear approximation of a hyperbolic model under small RR changes, derived similarly to the QT interval correction (QTc)38,39, Let’s call R Rr the reference RR interval associated to a reference heart beat and R R i the one to the i-th RR interval from one beat at the i-th segment, then As the Q T i − Q Tr difference, also dw(i ) is a measure of width change between the reference r and the current i-th mean T-waves from their respective observations time windows, then it is possible to extend previous relation in (11) to dw(i ) obtaining the HR related component By substituting (12) in (9) we obtain The value dw , c(i ) can be assumed to be non-zero mean, and uncorrelated to HR, that is: with dw , c(i ) zero mean and uncorrelated to HR. Then, dw(i ) becomes: where the parameters b, β and α , once jointly estimated (i.e. ˆ b , ˆ β and ˆ α ) can be used to derive ˆ dw , c(i ) as: Note that, ˆ β and ˆ α cannot be assessed from (15) with a directly least square fitting, since the DC component b in (15) largely affects the results. Rather, it is possible to jointly estimate ˆ b, ˆ β and ˆ α , and then use the results in (16). This estimate can be further approximated linearly for small RR changes. Denoting RR ( i ) =RRi−RRr , R R i can be expressed as R Ri=RRr+RR ( i ) and by replacing this in the right side of (12): Operating on the terms and under the assumption that RR ( i ) RRr , ( RR(i) RR r ) 1 and by using the Taylor’s series expansion, we have Substituting (18) in (13): where b, α , β and ( RRr) ( α−1 ) are constant values; then placing: the actual dw(i ) dependency with RR will be: (9) dw(i)=dw , c(i)+dw , HR(i) . (10) QT =β(RR)α. (11) Q Ti−QTr=β (RRi)α−(RRr)α . (12) d w,HR(i)=β (RRi)α−(RRr)α . (13) d w(i)=dw,c(i)+β (RRi)α−(RRr)α . (14) dw , c(i)=b+dw , c(i) , (15) dw(i)=b+dw,c(i)+β (RRi)α−(RRr)α , (16) ˆ d w,c(i)=dw(i)−ˆ β(RRi)ˆα−(RRr)ˆα . (17) dw,HR(i)=β (RRr+RR(i))α−(RRr)α . (18) (RRr+RR(i))α−(RRr)α≃αRR(i)(RRr)(α−1). (19) dw(i)≃b+dw,c(i)+αβRR(i)(RRr)(α−1), (20) αβ(RRr)(α−1)=c , (21) dw(i)≃b+dw , c(i)+cRR(i) .
; Vol:.(1234567890) Ƥ | (2021) 11:3883 | ǣȀȀǤȀͷͶǤͷͶ;ȀͺͷͻͿ;ǦͶͷǦ;ͿͻǦͻ www.nature.com/scientificreports/ From the geometrical point of view, b and c can be estimated as the zero-crossing and the slope, respectively, of the least-squares line fit to the dw(i ) values (in a RR ( i ) vs. dw(i ) graph). Then, the dw(i ) component that does not dependent on the RR, meaning it is assumed not correlated, can assessed as: where ˆc is the estimated slope from the Holter recording, see Fig.3, ˆ dw , c(i ) is then the corrected estimated of dw , c(i ) , with R R i and R Rr the mean RR interval from the i-th studied segment and the reference windows respectively and ˆc parameter is estimated for every patient during the time course of the Holter recording. When the linear approximation presented above cannot be assumed, ˆ b , ˆ β , and ˆ α can be jointly estimated from the model in (15), and use the (16) as the corrected estimate. An example of the estimated ˆ dw , c(i ) is given in Fig.5a where both ˆ dw , c(i ) and dw(i ) where displayed. Notice how the proposed correction formula removed the HR-dependency, for example around h = 12. Table3 provides an overview of the morphology markers studied in this work. [K+]Ǥ The proposed biomarkers have been compared with the relative variations in [ K+ ] (denoted as [K+](h ) ) with respect to a reference [ K+ ] that was taken at the end of the HD: (22) ˆ d w , c(i)=dw(i)−ˆcRR(i)=dw(i)−ˆc(RRi−RRr) . Table 3. T-wave morphology markers for [ K+ ] monitoring. *du w correspond to the marker denoted as d w in17, while here. d w is reserved for the newly introduced signed version. Markers Description Original markers from17 d u w *Time-domain changes between The reference and the i-th MWTW (ms). dN L w Nonlinear component of the Time-domain changes between The reference and the i-th MWTW (ms) da Relative amplitude changes between The reference and the i-th MWTW (%) dN L a Relative nonlinear amplitude changes After normalising the reference And the i-th MWTW (%) Specifically proposed in this work dw Signed version of the Previously proposed d u w * (ms) d w , cHeart rate corrected version of dw (ms) (b)(a) Figure3. Scatterplot showing the values of both d w panel (a) and ˆ dw , c panel(b) with respect to RR for a given patient in PCA approach. Spearman’s correlation coefficients ( ρ ) and p-values for both d w and ˆ dw , c are shown on top of each panel, while the least-square fitting regression lines are plotted in red.
Ϳ Vol.:(0123456789) Ƥ | (2021) 11:3883 | ǣȀȀǤȀͷͶǤͷͶ;ȀͺͷͻͿ;ǦͶͷǦ;ͿͻǦͻ www.nature.com/scientificreports/ being [ K+]h the concentration at the h-th hour during the HD and [ K+]r the concentration at the end of the treatment. An example of the [K+](h ) evolution is shown in Fig.5a (purple diamonds). Ǥ Results are presented as median and interquartile range (IQR). Spearman rank correlation coefficient ( ρ ) and Pearson correlation (r) were used for correlation analysis between [K+ ] and the proposed biomarker, giving information about both the monotonic relation and the strength of the association between the time warping based biomarkers and [ K+ ] changes and then providing a more complete characterisation. The average duration of the ECG recordings was 44 h mainly due to electrode detachment or early battery exhaustion. For this reason, correlation coefficients were computed using the first five values of [K+](h ) throughout the HD and the warping markers evaluated at the corresponding i-th segment points ( h=(i−1) /2 where i=1, 3, 5, 7, 9 or i=1, 3, 5, 7, 8 depending on the HD duration). All statistical analyses were performed using MATLAB version R2018b. In this study, ECG signals and [ K+ ] from 29 ESRD-HD patients were investigated. An example of d w and ˆ dw , c time evolution for a particular patient, in PCA approach, was provided in Fig.3. RR was represented on the x-axis in both panels, while d w and ˆ dw , c were shown on the y-axis in panel (a) and panel (b), respectively. The least-square fitting line (red line) was depicted in both panels. Spearman’s correlation coefficients ( ρ ) and p-values were also showed in each panel. High and significant correlation ( ρ=−0.9 0 and p-value < 0.00 1 ) was found between RR and d w . However, after correcting for the HR-dependency, ρ=0.0 3 and p-value =0.7 6 . Correlation between [ K+ ] and mean HR expressed in beats per minute (bpm) have also been computed and the results are presented in Table2, with a Spearman’s correlation coefficient median (IQR) values of 0.10 (1.35), and a median p-value of p=0.33. These values were 0.09 (1.45), p = 0.22 for Pearson’s correlation coefficient. Table4 shows the intra-patient Spearman’s ( ρ ) and Pearson’s (r) correlation coefficients computed between the relative variations in [ K+ ] (denoted as [K+ ] ) with respect to a reference [ K+ ] that was taken at the end of the HD and the time-warping parameters. In both single-lead and PCA approaches, the highest median Spearman’s and Pearson’s correlation coefficients were found for du w , d w and dw , c being ρ≥0.8 2 and r ≥0.8 6 for single-lead analysis and ρ≥0.8 2 and r ≥0.8 9 in PCA. Boxplots in Fig.4 show the distributions of [K+ ] and the proposed PCA-based time-warping descriptors during HD. Figure5b shows the average time evolution of PCA-based du w , d w , ˆ dw , c and dNL w in the studied population along the monitoring period, while the evolution of d a and dNL a is shown Fig.5c. Repolarization abnormalities play a fundamental role in the genesis of arrhythmic events and the risk increases in patients at ESRD with imbalance in [ K+ ] 41. In this work, two previously reported potassium estimators, Tw 14 and T S / A15, four warping-based ECG-derived biomarkers for [ K+ ] monitoring proposed in Ramírez etal.17, du w , d a , dNL w , dNL a , and the here proposed modified versions d w and dw , c were tested as bloodless indices for [ K+ ] variations in ESRD-HD patients computed from standard leads as well as in a PCA-derived lead. The most promising results in terms of correlation were obtained for markers du w , d w , and dw , c , leading to the highest median intra-patient ρ≥0.8 2 and r ≥0.8 7 in single-lead and ρ≥0.8 2 and r ≥0.8 9 in PCA lead respectively, evidencing high monotonic and linear association with [ K+ ] and making them a promising non-invasive indices for blood [ K+ ] monitoring. The signed biomarker d w followed a similar time-course as the unsigned du w during the whole monitoring period, showing a similar distribution in Fig.5a, as a result of the fact that the sign computed as in (5) is positive in roughly all the patients. That can be explained by the fact that the T-wave morphology in hyperkalemia is usually more peaked and shorter in time than a T-wave from regular [ K+ ] concentrations, as happens at the end of HD, where the reference has been taken42,43. Therefore, all the other MWTWs needed to be shrunk in amplitude and widened in time duration during the warping procedure to fit the reference one, and this is given by a positive signed d w . However, other external factors, such as the potassium removal rates44 or the dialysate potassium level45,46, might also have played a role in altering ventricular repolarization activity. The warping algorithm is applied over the MWTWs computed from different observing windows with different HRs, as is evident in Fig.5a. Therefore, a corrected version of the d w , derived similarly to the QT correction formula38,39, was proposed since the HR influences this marker as pointed out in Ramírez etal.17, and can be observed in Fig.5a as an example around hours h = 9, 12 and 43. A large number of models have been proposed for the computation of QTc values independent of HR37–40. However, a previous study38 found that the linear regression model fits better than any other model to the relationship between QT and the RR intervals. Also, for small RR variations, in section “Heart-rate-corrected T-wave warping” it is shown that hyperbolic QT to RR dependency becomes linear. Therefore, we used a linear model to derive an HR-corrected index, dw , c . This approach was used to estimate the d w component strictly related to [ K+ ] removing its relation with HR as showed in Fig.3, where the HR-dependency, clearly visible in panel (a), was cancelled after the correction, panel (b). Comparing the results for du w , d w and ˆ dw , c , all of them have proved to be highly correlated with [ K+ ] variations. However, it is important to remember that du w (and so d w ) is biased by the HR effects as previously described17, while ˆ dw , c is no longer dependent on it, possibly being responsible for the lower IQR in the correlation, 0.25, as compared to 0.35 and 0.36 for d w and du w , respectively (see Table4, PCA column). It should also be noted that the small differences between the ρ and r computed for ˆ dw , c and d w could be due to the low HR variations (23) [K+] ( h ) = ( [K+] h −[K+]r )