ORIGINAL RESEARCH published: 09 July 2019 doi: 10.3389/fphys.2019.00841 Risk Stratification in Idiopathic Dilated Cardiomyopathy Patients Using Cardiovascular Coupling Analysis Javier Rodriguez 1,2 , Steffen Schulz 3*, Beatriz F. Giraldo 1,2,4 and Andreas Voss 3 1 Institute for Bioengineering of Catalonia, The Barcelona Institute of Science and Technology, Barcelona, Spain, 2 Automatic Control Department (ESAII), Barcelona East School of Engineering (EEBE), Universitat Politècnica de Catalunya, Barcelona, Spain, 3 Institute of Innovative Health Technologies, Ernst-Abbe-Hochschule Jena, Jena, Germany, 4 Centro de Investigación Biomédica en Red de Bioengenieria, Biomateriales y Nanomedicina, Madrid, Spain Cardiovascular diseases are one of the most common causes of death; however, the early detection of patients at high risk of sudden cardiac death (SCD) remains an issue. The aim of this study was to analyze the cardio-vascular couplings based on heart rate variability (HRV) and blood pressure variability (BPV) analyses in order to introduce new indices for noninvasive risk stratification in idiopathic dilated cardiomyopathy patients (IDC). High-resolution electrocardiogram (ECG) and continuous noninvasive blood pressure (BP) signals were recorded in 91 IDC patients and 49 healthy subjects Edited by: Zbigniew R. Struzik, (CON). The patients were stratified by their SCD risk as high risk (IDCHR ) when after two The University of Tokyo, Japan years the subject either died or suffered life-threatening complications, and as low risk Reviewed by: (IDCLR ) when the subject remained stable during this period. Values were extracted from Luca Faes, ECG and BP signals, the beat-to-beat interval, and systolic and diastolic blood pressure, University of Palermo, Italy Philip Thomas Clemson, and analyzed using the segmented Poincaré plot analysis (SPPA), the high-resolution joint Lancaster University, United Kingdom symbolic dynamics (HRJSD) and the normalized short time partial directed coherence *Correspondence: methods. Support vector machine (SVM) models were built to classify these patients Steffen Schulz Steffen.Schulz@eah-jena.de according to SCD risk. IDCHR patients presented lowered HRV and increased BPV compared to both IDCLR patients and the control subjects, suggesting a decrease in Specialty section: their vagal activity and a compensation of sympathetic activity. Both, the cardio -systolic This article was submitted to Computational Physiology and and -diastolic coupling strength was stronger in high-risk patients when comparing with Medicine, low-risk patients. The cardio-systolic coupling analysis revealed that the systolic influence a section of the journal on heart rate gets weaker as the risk increases. The SVM IDCLR vs. IDCHR model Frontiers in Physiology achieved 98.9% accuracy with an area under the curve (AUC) of 0.96. The IDC and the Received: 07 February 2019 Accepted: 19 June 2019 CON groups obtained 93.6% and 0.94 accuracy and AUC, respectively. To simulate a Published: 09 July 2019 circumstance in which the original status of the subject is unknown, a cascade model was Citation: built fusing the aforementioned models, and achieved 94.4% accuracy. In conclusion, Rodriguez J, Schulz S, Giraldo BF and Voss A (2019) Risk Stratification in this study introduced a novel method for SCD risk stratification for IDC patients based Idiopathic Dilated Cardiomyopathy on new indices from coupling analysis and non-linear HRV and BPV. We have uncovered Patients Using Cardiovascular some of the complex interactions within the autonomic regulation in this type of patient. Coupling Analysis. Front. Physiol. 10:841. Keywords: idiopathic dilated cardiomyopathy, heart rate variability, blood pressure variability, coupling analysis, doi: 10.3389/fphys.2019.00841 sudden cardiac death, risk stratification Frontiers in Physiology | www.frontiersin.org 1 July 2019 | Volume 10 | Article 841 Rodriguez et al. Cardiovascular Coupling and Risk Stratification INTRODUCTION phenomenon in low-risk ischemic heart disease patients, as a measure of autonomic function. It is capable of predicting According to the 2015 update of the heart disease and stroke SCD-related mortality by assessing the absence of this behavior statistics of the American Heart Association (Mozaffarian et al., (Schmidt et al., 1999). The dynamics of the cardiovascular 2015), ∼325,000 cases of sudden cardiac death (SCD) occurred system behave in a highly complex way through the interplays in the United States in that year, and it is the cause of 15– of different linear and non-linear subsystems (Voss et al., 2009). 20% of mortality worldwide (Saour et al., 2017). The implantable Changes in blood pressure are reflected in changes in heart-rate cardioverter defibrillator (ICD) is commonly recommended in regulation, and vice versa (Cohen and Taylor, 2002). patients that are at high risk of suffering SCD, and the risk of Several linear and non-linear time series analysis approaches SCD is halved when one is implanted, although the presence or have been developed for the quantitative analysis of the absence of an ICD implant has no significant influence over the cardiovascular system in bivariate ways. However, linear rate of death itself (Kober et al., 2016). approaches might be insufficient to quantify non-linear The implantation of an ICD is suggested in patients with structures and the complexity of physiological systems. an ejection fraction (EF) lower or equal to 35%. Currently, Therefore, approaches from non-linear time-series analysis seem there is no effective means to stratify SCD risk in patients with to be more suited to capture complex interactions between EF above the risk threshold, who constitutes at least 70% of time series, and are able to quantify direct interrelationships the patients who will suffer SCD (Chugh, 2017). Additionally, such as the nonlinear influence of blood pressure on heart the effectiveness of ICD therapy is time dependent, making a rate. These coupling approaches are used to quantify direct prediction of the duration of the treatment desirable for the and indirect relationships, as well as causal and non-causal purpose of optimizing costs. Therefore, there is still a need for relationships between time series, providing deeper insights into additional predictors to identify patients with idiopathic dilated alterations of the cardiovascular system and leading to improved cardiomyopathy (IDC) who have an increased risk of SCD knowledge of the interacting regulatory mechanisms under and who benefit from ICD implantation (Bristow et al., 2004; different physiological and pathophysiological conditions. These Duray et al., 2006). approaches represent promising tools for detecting multivariate Previous studies have hypothesized various strategies to assess information flows (Schulz et al., 2013a). Some studies are based SCD risk based on clinical tests and data, imaging techniques, on the analysis of these interactions through the application of and signal processing methods, among others. The corrected bivariate coupling methods. QT interval was tested in elderly subjects and was associated For example, the directional cardiovascular interactions with SCD risk (Panikkath et al., 2011); others analyzed T-wave on young healthy subjects were assessed by bivariate and inversions, wide QRS-T angle, and the left bundle branch block multivariate coupling measures, finding that bivariate measures and found it prolonged the ability of QRS to predict all-cause better quantify the information transferred between indices, mortality, including SCD (Aro et al., 2011, 2012a,b), but despite while trivariate better reflects the existence and delay of directed being useful for prediction, it is unable to predict individual interactions (Javorka et al., 2017). Information decomposition risk. Sympathetic dominance of the autonomic nervous system measurements, in terms of variance or entropy, were explored in conjunction with pro-arrhythmic processes increases the to assess information dynamics in cardiovascular networks probability of SCD in patients with ventricular fibrillation (Faes et al., 2017), analyzing the heart period, the systolic problems (Schwartz et al., 1992). blood pressure, and the respiratory activity. The authors Studies related to cardiac imaging have analyzed the concluded that these measures of information transfer and myocardial scar, including the peri-infarction border zone, to information modification are better assessed through entropy- stratify SCD risk. The quantification of the total myocardial based and variance-based methods. Other work explored the scar was explored, proving to be superior to the left ventricular polysomnographic recordings and the finger blood pressure ejection fraction (LVEF) in predicting an appropriate ICD measurements in healthy subjects in order to investigate the therapy in cardiomyopathy patients (Bertini et al., 2012). Other differences between the wake-sleep states in the heart period studies demonstrated that the image-based analysis of myocardial and the systolic blood pressure coupling (Silvani et al., 2008). scars can contribute to the decision to implant ICDs in SCD They found that at low frequencies there are differences between patients (Schmidt et al., 2007; Roes et al., 2009; Wu et al., 2012). these states in human subjects. Additionally, the complexity and Heart rate variability (HRV) has been studied as a causality of the interactions of cardiovascular variability series measurement of autonomic tone. Higher vagal tone activity was assessed through linear model-based and non-linear model is related to increased spontaneous variations in heart rate, free techniques (Porta et al., 2014), and deducing that model and multiple non-linear techniques have been applied to free methods provide additional insights compared to the simpler study it (Wolf et al., 1978; Task-Force-of-the-European- linear model based approaches. Society-of-Cardiology-the-North-American-Society-of-Pacing- The behavior of cardiovascular coupling differs based on Electrophysiology, 1996; Ebrahimzadeh et al., 2014; Wu et al., physiological conditions. Consequently, we hypothesize that the 2014; Fujita et al., 2016). Lower levels of indices related to relationships between the cardiac and vascular systems will be this variability have been associated with patients at SCD different between the IDC patients at a high and low risk of SCD. risk, regardless of their LVEF (La Rovere et al., 1998, 2003). Therefore, the aim of this study is to analyze the suitability of Another index explored is heart rate turbulence, a consistent cardiovascular couplings for risk stratification in these patients. Frontiers in Physiology | www.frontiersin.org 2 July 2019 | Volume 10 | Article 841 Rodriguez et al. Cardiovascular Coupling and Risk Stratification We propose the characterization of these interactions through for cardiac sudden death (IDCHR ). None of these patients died features extracted from ECG and blood pressure signals, to better from a non-cardiac disease. Additionally, 49 healthy subjects (30 describe the complex dynamics of the cardiovascular interaction male, 19 female; aged 46 ± 14 years) were used as control group and identify new indices of cardiovascular risk. (CON). Table 1 presents the baseline clinical information of the IDC patients. In order to characterize autonomous regulation of the DATABASE cardiovascular systems, the following time series from ECG (RR interval) and BP signals were extracted using in- The signals analyzed in this study are part of German Autonomic house software (programming environment Delphi 3 and Regulation Trial (ART) study, oriented to evaluate risk predictors MATLAB R R2011b): of SCD and to improve the risk stratification in IDC patients. - Time series of heart rate consisting of successive beat-to-beat Signals from 220 IDC patients were recorded in two hospitals: the intervals [BBI, tachogram, (ms)] from the ECG signal. Friedrich-Schiller University Hospital in Jena (44 patients) and - Time series consisting of the maximum successive end-systolic the Franz-Volhard Clinic in Berlin (176 patients). All participants blood pressure amplitude values over time in relation to the provided their written informed consent to a protocol approved previous R-peak [SBP, systogram, (mmHg)] from the BP signal. by the local ethics committee of the two hospitals. This study - Time series consisting of the minimum successive end- complies with the Declaration of Helsinki. The study was diastolic blood pressure amplitude values over time in relation approved by the local Ethics Committee (No. 0986-11/02). to the previous R-peak [DBP, diastogram, (mmHg)] from the In the acquisition protocol, the ECG signals (32 bit resolution, BP signal. 1,600 Hz sampling frequency) and blood pressure (32 bit, 500 Hz) were synchronously recorded for 30 min. The ECG All extracted time series (BBI, SBP, DBP) were filtered by recordings were performed with a Porti system (TMSi BV, applying an adaptive variance estimation algorithm (Wessel et al., Netherlands), and the blood pressure recordings with the 2000) to remove and interpolate seldom occurring ventricular Portapress NIBP monitor (TNO Biomedical Instrumentation, premature beats and artifacts (e.g., movement, electrode noise, Amsterdam, the Netherlands). All patients were recorded under and extraordinary peaks) to obtain normal-to-normal beat resting conditions in supine position. time series (NN). To obtain synchronized time series, BBI, The inclusion criteria for the IDC patients were LVEF <45% SBP, and DBP were resampled using a linear interpolation and/or fractional shortening <25%, and the left ventricle end- method (2 Hz). diastolic diameter (LVEDD) >117%. Additionally, the New York heart association index (NYHA) of each patient was included. The exclusion criteria were systemic hypertension at METHODS rest, coronary artery disease, congenital heart disease, pericardial diseases, valvular heart diseases, systemic disease known to The BBI, SBP, and DBP time series were extracted using cause dilated cardiomyopathy, chronic alcoholism, sustained algorithms based on zero crossings and different thresholds. ventricular tachycardia, atrial fibrillation, diabetes mellitus, renal These data were evaluated through several non-linear failure, and active permanent cardiac pacemaker as well as characterization techniques such as high resolution joint patients without sinus rhythm (Voss et al., 2010, 2012b). symbolic dynamics (HRSJD), segmented Poincaré plot Afterwards, an additional exclusion criteria was applied to analysis (SPPA) and normalized short-time partial directed discard patients with comorbidities and confounding factors coherence (NSTPDC). influencing the autonomic regulation system. 129 of the 220 patients were excluded due to the following reasons: 36 with paced rhythm, 34 patients suffered coronary artery TABLE 1 | Baseline clinical information of ART database IDC patients (median and disease, 32 presented a high (>5%) percentage of ectopic interquartile range). beats or artifacts, 19 with atrial fibrillation, four because IDCHR IDCLR p-value of technical problems, two suffered from hypertrophic non- N = 77 (59 ♂, 18 ♀) N = 14 (11 ♂, 3 ♀) obstructive cardiomyopathy, one patient with arrhythmogenic right ventricular cardiomyopathy, and another patient who was Follow-up duration 27 [17; 37] 30 [21; 38] n.s. clinically unstable due to acute decomposition (Voss et al., [months] 2012b). Finally, a total number of 91 IDC patients (21 female, 70 Age (years) 55 [50; 60] 56 [50; 63] n.s. male) were investigated in this work. LVEF (%) 29 [27; 37] 35 [27; 46] n.s. After a median follow-up period of 28 months (range: LVEDD (mm) 69 [61; 79] 61 [58; 68] n.s. 17–38 months), the patients were classified into two groups LVESD (mm) 60 [53; 69] 49 [44; 56] 0.0028 according to their SCD risk. The group of patients that remained NYHA 3 [2; 4] 2 [2; 3] 0.0024 in stable physical condition were considered as low risk for IDCHR and IDCLR , high risk and low risk group of patients with dilated cardiomyopathy; cardiovascular sudden death (IDCLR ). The remaining patients, N, number of patients within the groups (♂ = Male; ♀ = Female); LVEF, Left ventricular who either died of SCD or needed resuscitation because of a ejection fraction; LVEDD, Left ventricular end-diastolic diameter; LVESD, left ventricular life-threatening tachyarrhythmia, were categorized as high risk end-systolic diameter; NYHA, New York heart association index. Frontiers in Physiology | www.frontiersin.org 3 July 2019 | Volume 10 | Article 841 Rodriguez et al. Cardiovascular Coupling and Risk Stratification High Resolution Joint Symbolic Dynamics the transformation from X to S varies as is shown in Equation (3). The joint symbolic dynamics (JSD) method (Baumert et al., 2002) BBI − xBBI < −lBBI   is based on the analysis of dynamic processes by the means   0 : xn+1 n of symbols. Considering BBI, SBP, and DBP time series, X is  BBI − xBBI ≤ lBBI 1 : −lBBI ≤ xn+1   defined as a bivariate sample vector that contains two out of the SBBI n = n  BBI − xBBI > lBBI  three time series, for all their possible combinations (BBI-SBP, 2 : xn+1    n BBI-DBP, and SBP-DBP), as expressed in Equation (1), SBP − xSBP < −lSBP     0 : xn+1 n  SBP ≤ xSBP − xSBP ≤ lSBP T    X BBI_SBP = xnBBI , xnSBP  SSBP = 1 : −l  n n+1 n (3) T  X BBI_DBP = xnBBI , xnDBP  n = 0, 1, . . . , N with x ∈ R, (1)  SBP − xSBP > lSBP  2 : xn+1    SBP DBP T   n X SBP_DBP = xn , xn  DBP − xDBP < −lDBP     0 : xn+1 n  DBP ≤ xDBP − xDBP ≤ lDBP   being N the total number of samples. SDBP n = 1 : −l n+1 n In JSD, the increments between two successive values of    DBP − xDBP > lDBP 2 : xn+1   n the temporal series are coded as “1” and the decrements and equilibriums are coded as “0,” these increments and decrements are considered in relation to a threshold l. The vector X can be The new space is composed of combinations of 27 different transformed into the symbolic vector S using the rules set out in possible types of words (from 000 to 222) and a total of 729 Equation (2), were the threshold l is set equal to 0 (Baumert et al., indices. All the word types were grouped into eight pattern 2002; Giraldo et al., 2015). families, transforming the vector matrix W into a vector matrix family (Wf ). These indices were analyzed by their occurrence probabilities. Afterwards, these indices were grouped T  SBBI_SBP = SBBI , SSBP according to their family description (Table 2). These pattern  n n  BBI_DBP  BBI DBP T  families represent different interactions between the branches S = Sn , Sn n = 0, 1, . . . , N, with S ∈ 0, 1 DBP T    of the autonomic regulation system, leading to indices with SSBP_DBP = SSBP  n , S n statistically sufficient occurrence probabilities (Schulz et al., BBI − xBBI ≤ 0 (  0 : xn+1 n 2015b). Additionally, the Shannon entropy (Rundle et al., 2019) BBI Sn = BBI − xBBI > 0  (2) was calculated for all the proposed families to assess the 1 : xn+1 n complexity of the coupling. SBP − xSBP ≤ 0 (  0 : xn+1 n The thresholds applied were 5 ms and 1 mmHg on the SBP Sn = SBP − xSBP > 0  BBI and the blood pressure time series, respectively. These 1 : xn+1 n threshold values were successfully applied in earlier work (Schulz DBP − xDBP ≤ 0 (  0 : xn+1 n et al., 2013b). The threshold level using spontaneous baroreflex DBP Sn = DBP − xDBP > 0  sensitivity, in contrast to other thresholds, is the most suitable for 1 : xn+1 n highlighting different specific cardiovascular coupling patterns. Segmented Poincaré Plot Analysis Sequences of symbols are considered words of length k. The The Poincaré plot analysis (PPA) is used to quantify self- words are then arranged in a vector matrix W. Particularly, for similarity in processes by plotting the data into a higher k = 3 (Sn , Sn+1 , Sn+2 ) an 8 × 8 vector matrix can be derived, dimensional state space. Considering a time series X (n) = taking values from (000, 000)T to (111, 111)T . In order to obtain JSD indices that are more robust against noise, fluctuations and artifacts, the comparison threshold should TABLE 2 | Description of pattern families explored in the HRJSD method. be something other than 0. There are several advantages for Family Description choosing a non-zero threshold. For instance, a state will be generated that will help to distinguish between small and E0 No variation of 3 successive “0” symbols (“000”) large changes in the system’s variability. It is also possible to E1 No variation of 3 successive “1” symbols (“111”) differentiate between decrements and equilibrium because both E2 No variation of 3 successive “2” symbols (“222”) states are no longer coded with the same symbols. Lastly, the LU1 Low increasing behavior (“122,” “022,” “112,” “221,” “220,” “211,” number of word types including “0” will not be the most “121,” “212”) prevalent within the W matrix (Schulz et al., 2013b, 2015b). LD1 Low decreasing behavior (“011,” “001,” “002,” “110,” “100,” The high-resolution joint symbolic dynamics (HRJSD) method “200,” “010,” “101”) implements a three symbol JSD after setting a threshold l. The LA1 Fast alternant behavior (“020,” “202”) increment states are coded as “2,” the decrements are coded as “1” P Alternant peak-like behavior (“120,” “201,” “210”) and the equilibrium states are coded as “0.” With this technique, V Alternant valley-like behavior (“021,” “102,” “012”) Frontiers in Physiology | www.frontiersin.org 4 July 2019 | Volume 10 | Article 841 Rodriguez et al. Cardiovascular Coupling and Risk Stratification x1 , x2 , x3 , . . . , xn , the Poincaré plot is obtained by plotting X (n) to yield information about the partial correlative short-time vs. X (n + 1). The typical representation is an elongated scatter of interaction properties (Baccala and Sameshima, 2001). plots through the line of identity with all the points whose values The NSTPDC is able to manage non-stationary signals by are near the mean placed toward the center. Long- and short- evaluating their dynamic coupling changes and detecting their term variability indices can be by fitting an ellipse to the shape of level and directions in multivariate and complex dynamic the plot and measuring the dispersion along the minor (SD1) and systems (Adochiei et al., 2013; Schulz et al., 2015a). major (SD2) axis of the ellipse, corresponding to their standard Normalized short-time partial directed coherence based on deviation (Rodriguez et al., 2017). an m-dimensional multichannel auto-regressive model (MAR) Although PPA is a non-linear characterization method, process with model order p to determine Granger causality in the indices SD1 and SD2 can be correlated with the linear behavior frequency domain. For the selection of the optimal model order of the system, hence making them a suboptimal way to explore popt of the AR(p) model and for the estimation of its coefficients, information about the non-linear part of the process (Seeck et al., the stepwise least squares algorithm (Neumaier and Schneider, 2011; Voss et al., 2012a). 2001) and the Schwarz’s Bayesian Criterion (SBC) were applied The Segmented Poincaré plot analysis (SPPA) is an enhanced (Schneider and Neumaier, 2001). NSTPDC is based on the time- pseudo-phase space quantification method that yields indices variant partial directed coherence approach (tvPDC, πxy (f , n)) that also represent the non-linear information of the system. In providing information about the partial correlative short-time SPPA the SD1 and SD2 indices are calculated similarly to PPA. interaction properties of non-stationary signals, with f as the Then, the scatter points are rotated α degrees around the main frequency and n the number of windows (Milde et al., 2011). focus of the plot as defined in Equation (4), To quantify the coupling direction between two time series, X and Y (e.g., BBI and SBP: with xBBI and ySBP ) with the covariate z  ′        (e.g., DBP with zDBP ), a coupling factor (CF) was introduced. CF Xn Xn cos α − sinα 0 Xn − Xn was obtained by dividing the mean value πxy (f , n) by the mean ′ X n+1  = Xn+1  +  sin α cos α 0  × Xn+1 − Xn+1  . ′ value of πyx (f , n), defined as z z 0 0 1 z (4) 1 P πxBBI ySBP f , n  1X n  CF = 1  , a= πxBBI ySBP f , n , n P πySBP xBBI f , n The x and y axis correspond to X (t) and X (t + 1) values, Xn and n Xn+1 are the mean values of the original and shifted X time series 1X  b= πySBP xBBI f , n . (6) respectively, and z is the axis of the rotation (Voss et al., 2010, n 2012a; Seeck et al., 2011). In this case we define α = 45 degrees in order to simplify the estimating procedure of the SD1/SD2 These results were normalized to become a specific set of values adapted probability. Afterwards, a 12 × 12 rectangular grid is leading to the (normalized) factor NF representing the coupling drawn for the plot. The size of the rectangles (height, width) is direction, given by adapted based on the SD1 (row) and SD2 (column) values.    Finally, for each rectangle in position (i, j) the single a   2, if max = a & > 5  b probability (ρij ) is calculated considering the number of points NF = 1, if (max = a & 2 < a ≤ 5) and contained by the total number of points in the series. Afterwards,  b  0, if (max = a & 0 ≤ a ≤ 2)  the probabilities of each row (ρri ) and column (ρcj ) are calculated b as the sum of their single probabilities, as shown in Equation (5):    b  −2, if max = b & a > 5    12 NF = −1, if max = b & 2 < b ≤ 5 . (7) X a (5)  ρri = ρij  0, if (max = b & 0 ≤ ab ≤ 2)  j= 1 12 X A normalization procedure was applied to CF leading to the ρcj = ρij . normalized factor (NF). NF determinates the direction of the i=1 causal connections between the investigated time series (xBBI and ySBP ) as a function of frequency f. NF takes the following Normalized Short-Time Partial values: NF = {−2, −1, 0, 1, 2}. Strong unidirectional coupling Directed Coherence is indicated if NF is −2 or 2 (where −2 denotes ySBP as the The directed coherence method (DC) describes how and whether driver), bidirectional coupling if NF = −1 or 1 (−1 denotes ySBP two complex physiological signals are functionally connected as the driver), and an equal influence in both directions and/or (Faes et al., 2012). The DC method studies the relative structural no coupling if NF = 0 in respect to coupling strengths (if both relationships between the systems by breaking down their area indices reveal equal values that are larger than zero an equal interactions into feedback and feedforward aspects. influence in both directions is present, on the other hand, if both The partial directed coherence method (PDC) determines the area indices reveal equal values but are zero there is no coupling). either direct or indirect causality between the systems analyzed. For determining the coupling strength between two time The PDC is limited to work on stationary signals and is unable series, e.g., xBBI and ySBP with covariate (zDBP ), the areas Frontiers in Physiology | www.frontiersin.org 5 July 2019 | Volume 10 | Article 841 Rodriguez et al. Cardiovascular Coupling and Risk Stratification (ABBI→SBP(DBP) , ASBP→BBI(DBP) , [a.u.]) generated in space by the entire recording) by the total number of NN intervals CF were estimated in each window within the frequency band (BBI_PNN50, SBP_PNN50, DBP_PNN50). f = 0–2 Hz, and afterwards, averaged. ABBI→SBP(DBP) and - The square root of the mean squared differences of BBI, SBP, ASBP→BBI(DBP) ranges between 0 to 1 [0, 1]. Hereby, 1 indicates and DBP time series (BBI_rmssd, SBP_rmssd, DBP_rmssd). that all causal influence originating from time series X are - The power of the low frequency components (0.04–0.15 Hz) of directed toward (arrows:→ ) time series Y. BBI, SBP, and DBP time series (BBI_LF, SBP_LF, DBP_LF). In order to take advantage of the aspect of stationarity and - The power of the high frequency components (0.15–0.4 Hz) of scale-invariance for NSTPDC analyses, a normalization (zero BBI, SBP, and DBP time series (BBI_HF, SBP_HF, DBP_HF). mean and unit variance) of the time series BBI, SBP, and DBP - The ratio between the low and high frequency power was performed (Schulz et al., 2015a). Therefore, each sample i of components of BBI, SBP, and DBP time series (BBI_LF/HF, the BBI-,  SBP-, or DBP- time series X = {xi , i = 1, . . . N} and SBP_LF/HF, DBP_LF/HF). Y = yi , i = 1, . . . N with N as the maximal number of samples i (temporal index) was first normalized by subtracting the mean Feature Extraction of x and, then divided by the standard deviation (std) of X or A total number of 621 indices were extracted to analyze the Y, respectively. The normalized time series xnorm and ynorm (zero interaction between the coupling across all signals (ECG and mean and unit variance) were thus obtained as seen in (8) blood pressure signals). The “cd” label represents the cardiac and x (i) − x y (i) − y diastolic coupling, “cs” label indicates the cardiac and systolic xnorm (i) = and ynorm (i) =  . (8) coupling, and “ds” label the diastolic and systolic coupling. std (x) std y The distribution of these indices is as follows: 264 indices from Dual Sequence Method the HRJSD, 216 from the SPPA, 97 from the JSD, 12 from the PPA, A widely spread method to investigate the spontaneous 9 from the NSTPDC, 21 from the standard HRV and BPV indices, baroreflex sensitivity (BRS) is the sequence method (Bertinieri and 2 from the DSMs. Summary descriptions of the indices are et al., 1988). The BRS was obtained scanning the SBP and BBI shown in Table 3. time series for sequences of three or more successive heart beats in which a progressive increase (or decrease) in SBP was followed (with a one-beat delay). The dual sequence method (DSM) was STATISTICAL ANALYSIS developed to improve the analysis of the baroreflex sensitivity In order to reduce dataset dimensionality, Mann Whitney non- (Malberg et al., 1999, 2002). A fluctuation is defined when there parametric statistical test was used to determine the statistical are variations greater than 1 mmHg (increasing or decreasing) significance of the indices obtained in the characterization in SBP and greater than 5 ms in BBI values. The highest slope of every sequence was taken for the linear regression. The slopes of the regression lines between the SBP and BBI sequences were TABLE 3 | Coupling indices extracted from high-resolution joint symbolic taken as an index for local BRS [ms/mmHg] and calculated dynamics, segmented Poincaré plot analysis, and normalized short-time partial in each recording. The DSM is based on standard sequence directed coherence. methods, the enhancement lies in the analysis of two different Index Description kinds of BBI response: bradycardic (an increase in SBP causes an increase in BBI) and tachycardic fluctuations (a decrease in HRJSDxy_Fx-Fy Probability of occurrence of the Fx and Fy word SBP causes a decrease in BBI), whereas only the bradycardic families from the x–y coupling fluctuations represented the classical spontaneous baroreflex HRJSDShxy Shannon entropy of all the word families from the sensitivity (bslope). The analysis of the tachycardic fluctuations x-y coupling (tslope) provides additional information about autonomous HRJSDxy-Fx Summation of the occurrences of Fx on all the families from the x-y coupling cardiovascular regulation (Parati et al., 2000). SPPAxy_Row_n-m Probability of occurrence of the Row n-m from the Heart Rate and Blood Pressure Variability x-y coupling SPPAxy_Column_n-m Probability of occurrence of the Column n-m from Standard Indices the x-y coupling The HRV and the blood pressure variability (BPV) were JSDxy-n Probability of occurrence of the word n from the x-y quantified using the standard indices from time and frequency coupling domain. The following indices were calculated from the BBI, PPAxy_SD1 Short-time standard deviation from the x-y coupling SBP, and DBP time series, with NN being the normal-to- PPAxy_SD2 Long-time standard deviation from the x-y coupling normal intervals: PPAxy_SD1/SD2 Short and long deviations ratio from the x-y coupling - The mean of BBI, SBP, and DBP time series (BBI_meanNN, NSTPDCxy_NF Normalized coupling factor from the x-y coupling SBP_meanNN and DBP_meanNN, respectively). NSTPDCxy_ Ax→ y Coupling strength of x-y from x to y - The standard deviation of BBI, SBP, and DBP time series NSTPDCxy_ Ay→ x Coupling strength of x-y from y to x (BBI_sdNN, SBP_sdNN, DBP_sdNN). x and y represent the couplings cd, cs, and ds; F represents the word families: E0, E1, - The proportion derived by dividing NN50 (number of pairs E2, LU1, LD1, LA1, P and V; n and m represent the number of row and column in the of adjacent NN intervals differing by more than 50 ms in PPA based indices, respectively. Frontiers in Physiology | www.frontiersin.org 6 July 2019 | Volume 10 | Article 841 Rodriguez et al. Cardiovascular Coupling and Risk Stratification process. The results were analyzed for different levels of The classification problem is solved by maximizing the margin significance, including the Bonferroni criterion, considering: while minimizing the training error. Using the Lagrange multipliers method, a dual formulation can be obtained (Cortes and Vapnik, 1995) (Equation 13), n.s. ρ ≥ 0.01 ∗ ρ ≤ 0.01 significant 1 X minP w, b = kwm zk2 + C    ∗∗ ρ ≤ 0.001 highly significant K1 yi f (xi ) (13) ∗∗∗ ρ ≤ 0.0000167 Bonferroni criterion 2 i n = 621 indices where C is a penalty parameter. Besides the scale of C having Additionally, a correlation analysis was performed on the no direct meaning, as its value increases, the penalty assigned statistically significant indices. The ones with high correlation to errors is stronger, narrowing the decision boundary (Ben-Hur (ρ ≥ 0.7) and relative lower significance were discarded. et al., 2008). The leave-one-out cross-validation procedure was used to Each feature was scaled and normalized (zero mean and unit validate the results. The classification results are presented in variance) in order to avoid scaling biases. For each iteration of terms of accuracy (Acc), sensitivity (Sn), specificity (Sp), and area features, the model was built by optimizing the value of C for each under the curve (AUC). of the kernels considered, by iterating different values of σ and d. The indices that showed statistical differences and low correlation CLASSIFICATION TECHNIQUES were used in pairs to build several SVM models. The accuracy of each model was then calculated and the one with the higher value The objective of the support vector machines (SVM) method is was chosen as optimal for each type of kernel. to find a higher dimensional space where a classification problem can be more easily solved than in the original space. The vectors RESULTS that defines the hyperplane are called support vectors. This technique allows separating groups of data that are not originally The calculated indices were used to analyze the cardiovascular possible using linear classifiers (Giraldo et al., 2015). coupling in 91 IDC patients and 49 healthy subjects. Four Being  X = {x1 , . . . , xL } , x ∈ R a given set of data vectors and different comparisons were performed: Y = y1 , . . . , yL their corresponding labels, the SVM function, • The high risk IDC patients (IDCHR ) vs. the low risk IDC defined as a linear discriminant function, known as hyperplane, patients(IDCLR ). is given can be defined by Equation (9), • The IDC patients vs. the CON subjects. XL  • The high risk IDC patients (IDCHR ) vs. the CON subjects. f (x) = wz + b = αi yi K xi yi + b, (9) • The low risk IDC patients (IDCLR ) vs. the CON subjects. i where w is the normal vector to the hyperplane. The function IDCHR Patients vs. IDCLR Patients K(xi yi ) is the Kernel function that will shape the hyperplane and Compared With CON Subjects αi and b define the efficiency of the classifier on the optimal We obtained statistically significant differences in the symbolic values. In this study we evaluated the Gaussian, Laplace, and dynamic analysis, for both the cardio-diastolic (BBI–DBP) and ANOVA kernels. diastolic-systolic (DBP-SBP) couplings. Figures 1, 2 present an The Gaussian kernel is useful to model radially distributed example of the three-dimensional plots of the word distribution data, and is defined in Equation (10), where σ is a density matrix of the couplings, using the HRJSD method from penalization term, (Figures 1A,C, 2A,C) IDCLR and (Figures 1B,D, 2B,D) IDCHR patients, respectively.  kx−yk2  Figure 3 shows an example of the Poincaré plot method − 2 applied to a patient for each analyzed group, considering the  2σ K x, y = e . (10) systogram from (Figure 3A) a CON subject, (Figure 3B) a IDCLR Laplace kernel is a less σ influenced version of the Gaussian patient, and (Figure 3C) a IDCHR patient. kernel, given by Equation (11) Figure 4 represents the averaged NSTPDC applied on the BBI, SBP, and DBP time series couplings, for (Figure 4A) the CON,  kx−yk2  (Figure 4B) IDCLR , and (Figure 4C) IDCHR groups. −  K x, y = e 2σ . (11) Figure 5 presents the relationships between all three analyzed groups, where the arrows represent the coupling direction and the arrow thickness indicating the coupling strength. In addition, The ANOVA kernel works well on multidimensional support the level of statistical significance has been represented (p ≤ 0.01). vector regression models (Stitson et al., 1999) and is defined by When comparing the IDCHR and IDCLR groups, 96 indices Equation (12), where σ and d are the optimization indices of presented statistically significant differences, corresponding to: its function, 76 indices from the HRJSD, 12 from the SPPA, and 7 from the  Xn k −yk )2 ) d NSTPDC. After the correlation analysis, a total of 36 statistically K x, y = e(−σ (x . (12) significant indices were chosen for the classification process. k=1 Frontiers in Physiology | www.frontiersin.org 7 July 2019 | Volume 10 | Article 841 Rodriguez et al. Cardiovascular Coupling and Risk Stratification FIGURE 1 | Three-dimensional plots of the word distribution density matrix using the HRJSD method (single word probabilities, word families), from a patient of (A,C) IDCLR and (B,D) IDCHR , respectively, for cardio-diastolic coupling. Some of the most relevant indices, expressed in a mean value and HRV and BPV Standard Indices and 95% confidence interval, are shown in Table 4. DSM Results When the HRV and BPV standard indices and the DSM IDC Patients vs. CON Subjects were evaluated in the IDCHR vs. IDCLR comparison, no When comparing IDC patients vs. CON subjects, 261 indices statistically significant indices were found. In the remaining presented statistically significant differences between the two comparisons, seven indices presented statistical significances: groups: 169 from the HRJSD, 76 from the SPPA, and 13 two from the sequence analysis and five from the HRV and from NSTPDC. A total of 82 indices remained after discarding BPV standard indices. A summary of these results is shown the highly correlated indices. Table 5 shows the most relevant in Table 8. of them. Classification Results IDCHR and IDCLR Patients vs. After the SVM classification step, the PPAs_SD1/SD2 and CON Subjects HRJSDds_LU1-P indices were deemed the optimal choices for When comparing IDCHR patients vs. CON subjects and IDCLR the Laplace kernel SVM model, achieving an accuracy of 98.9% patients vs. CON subjects, the differences were found in 182 and and an AUC of 0.96 for the IDCHR vs. IDCLR comparison. 247 indices, respectively. After the correlation analysis, 61 and The HRJSDds-E0d and NSTPDCcs_Area_c-s allowed us to 85 remaining indices were chosen for the classification step. A classify IDC patients from the CON group with an accuracy summary of the most relevant indices is shown in Tables 6, 7. of 93.6 % and an AUC of 0.94 using the Laplace kernel. Frontiers in Physiology | www.frontiersin.org 8 July 2019 | Volume 10 | Article 841 Rodriguez et al. Cardiovascular Coupling and Risk Stratification FIGURE 2 | Three-dimensional plots of the word distribution density matrix using the HRJSD method (single word probabilities, word families), from a patient of (A,C) IDCLR and (B,D) IDCHR , respectively, for diastolic-systolic coupling. A B C FIGURE 3 | Systolic blood pressure Poincaré plot analysis results from (A) a CON subject, (B) a IDCLR patient, and (C) a IDCHR patient. Frontiers in Physiology | www.frontiersin.org 9 July 2019 | Volume 10 | Article 841 Rodriguez et al. Cardiovascular Coupling and Risk Stratification FIGURE 4 | Averaged NSTPDC plots for cardiovascular coupling analyses for (A) the CON, (B) the IDCLR , and (C) the IDCHR group. Arrows indicating the causal coupling direction from one time series to another time series, e.g., SYS←BBI, indicating the causal link from BBI to SYS. Coupling strength ranges from blue (0, no coupling) to red (1, maximum coupling) where BBI are beat-to-beat intervals, and SYS are successive end-systolic blood pressure amplitude values over time. FIGURE 5 | Graphical representation of the entire cardiovascular coupling structure (coupling strengths and directions) among the cardiac (BBI), systolic blood pressure (SYS), and diastolic blood pressure (DIA) systems, when comparing (A) CON vs. IDCLR , (B) CON vs. IDCHR , and (C) IDCHR vs. IDCLR . The arrows direction indicates the causal coupling direction and the thickness the coupling strength. *Indicates that the NSTPDC index associated to the coupling represented by the arrow is statistically significant when *p ≤ 0.01; **p ≤ 0.001; ***p ≤ 0.0000167. Meanwhile, the HRJSDds_LD1-V and the NSTPDCcs_Area_c- The cascade model achieved 94.4% accuracy. From a One s were able to discriminate between IDCHR patients and the vs. All approach, the results for the CON group were 86.1% CON group with an accuracy of 96.9% and AUC of 0.95 sensitivity and 92.8% specificity, from the perspective of the HR applying the Gaussian kernel. Finally, the HRJSDds-E1d and group, the model yielded a sensitivity of 93.2 and specificity of the NSTPDCcs_Area_c-s were found to be the best indices to 94.5%, and for the LR group, the model achieved 98.8% sensitivity classify IDCLR patients from the CON group, obtaining 89.6 and 89.0% specificity. % accuracy and a 0.85 AUC with the Laplace kernel. The classification plots and the results are shown in Figure 6 and Table 9, respectively. DISCUSSION Cascaded Risk Stratification The aim of this study was to find indices capable of stratifying In order to consider the typical clinical case in which the original SCD risk in IDC patients. In order to achieve this, we derived condition of the subject is unknown, a cascade model was indices extracted from BBI (ECG) and systolic and diastolic developed to classify this new subject using the label SDC risk, blood pressure (BP) temporal series using linear and non-linear comparing the IDC vs. CON and IDCHR vs. IDCLR models univariate and bivariate (coupling analysis) techniques. The (Figure 7). The general idea of this model is to classify the subject indices extracted from these techniques revealed patterns that in either CON, LR, or HR without any prior labeling. The first behave differently in patients at high risk of SCD. These indices step is to decide if a subject is an IDC patient or not using the IDC (mainly from coupling analyses) were used to train several SVM vs. CON model. Afterwards, those classified as IDC patients are models in order to classify the subjects based on different levels analyzed based on the LR vs. HR model to predict their risk level. of SCD risk. Frontiers in Physiology | www.frontiersin.org 10 July 2019 | Volume 10 | Article 841 Rodriguez et al. Cardiovascular Coupling and Risk Stratification TABLE 4 | Mean value and 95% confidence interval of the most significant indices TABLE 6 | Mean value and 95% confidence interval of the most significant indices when comparing the IDCHR and IDCLR patients. comparing the IDCHR patients and CON subjects. Index IDCHR (N = 14) IDCLR (N = 77) p-value Index IDCHR (N = 14) CON (N = 49) p-value PPAs_SD1/SD2 0.26 [0.22; 0.30] 0.20 [0.18; 0.22] * HRJSDds_LD1-V 0.05 [0.040; 0.066] 0.02 [0.020; 0.027] *** PPAd_SD1/SD2 0.27 [0.25; 0.28] 0.19 [0.17; 0.20] ** NSTPDCcs_Area_c->s 0.29 [0.24; 0.35] 0.15 [0.13; 0.17] *** PPAd_SD1 1.21 [1.15; 1.26] 0.92 [0.85; 0.98] ** NSTPDCcd_Area_c->d 0.31 [0.27; 0.36] 0.22 [0.20; 0.24] ** SPPAcd_Column_5 13.97 [13.54; 14.39] 12.09 [12.01; 13.02] * SPPAds_Column_2-6 27.23 [25.73; 28.72] 22.44 [21.26; 23.63] ** SPPAcd_Column_8 13.77 [13.53; 14.03] 12.71 [12.35; 13.06] n.s. HRJSDcd_LD1-E2 0.003 [0.001; 0.006] 0.012 [0.010; 0.015] ** HRJSDcd_E1-LD1 0.03 [0.024; 0.035] 0.09 [0.081; 0.091] ** **ρ ≤ 0.001; ***ρ ≤ 0.0000167; IDCHR , idiopathic cardiomyopathic patients at high HRJSDcd_ LD1-V 0.056 [0.039; 0.073] 0.025 [0.021; 0.028] ** risk of sudden cardiac death; CON, control group; SPPA, segmented Poincaré plot HRJSDds_P-LU1 0.028 [0.022; 0.034] 0.016 [0.013; 0.019] ** analysis; HRJSD, high resolution joint symbolic dynamic; NSTPDC, normalized short-time HRJSDds_V-LD1 0.050 [0.040; 0.060] 0.026 [0.022; 0.029] ** partial directed coherence; cd, cardio-diastolic coupling; cs, cardio-systolic coupling; ds, diastolic-systolic coupling. NSTPDCcd_NF −0.15 [−0.44; 0.14] −0.71 [−0.85; −0.56] * NSTPDCcd_Area_c-d. 0.26 [0.23; 0.29] 0.20 [0.18; 0.21] ** *ρ ≤ 0.01; **ρ ≤ 0.001; n.s. ρ ≥ 0.01; IDCHR , idiopathic cardiomyopathic patients at high TABLE 7 | Mean value and 95% confidence interval of the most significant indices risk of IDCLR , idiopathic cardiomyopathic patients at low risk of SCD; PPA, Poincaré plot comparing the IDCLR patients and CON subjects. analysis; SPPA, segmented Poincaré plot analysis; HRJSD, high resolution joint symbolic dynamic; NSTPDC, normalized short-time partial directed coherence; cd, cardio-diastolic Index IDCLR (N = 77) CON (N = 49) p-value coupling; ds, diastolic-systolic coupling. HRJSDcd- 0.33 [0.29; 0.37] 0.21 [0.17; 0.26] ** E1d TABLE 5 | Mean value and 95% confidence interval of the most significant indices comparing IDC patients and CON subjects. NSTPDCcs_ 0.25 [0.22; 0.27] 0.15 [0.13; 0.17] *** Area_c-s Index IDC (N = 91) CON (N = 49) p-value SPPAcd_ 5.40 [4.95; 5.85] 7.19 [6.56; 7.82] ** Column1–8 HRJSDds- 0.018 [0.014; 0.021] 0.0064 [0.0048; 0.0071] *** HRJSDcd_ 0.0004 [0.0002; 0.0005] 0.0014 [0.0010; 0.0019] *** E0d E0-LA1 NSTPDCcs_ 0.25 [0.23; 0.28] 0.15 [0.13; 0.17] *** HRJSDds_ 0.070 [0.064; 0.075] 0.093 [0.086; 0.101] *** Area_c-s LD1-LD1 PPAc_SD1 14.18 [10.44; 17.93] 23.27 [20.70; 25.84] ** **ρ ≤ 0.001; ***ρ ≤ 0.0000167; IDCLR , idiopathic cardiomyopathic patients at low PPAc_ 0.29 [0.23; 0.34] 0.36 [0.33; 0.38] ** risk of sudden cardiac death; CON, control group; SPPA, segmented Poincaré plot SD1/SD2 analysis; HRJSD, high resolution joint symbolic dynamic; NSTPDC, normalized short-time HRJSDcd- 0.02 [0.18; 0.21] 0.05 [0.045; 0.054] * partial directed coherence; cd, cardio-diastolic coupling; cs, cardio-systolic coupling; ds, E2d diastolic-systolic coupling. HRJSDcd- 0.23 [0.22;0.24] 0.28 [0.27; 0.29] ** LD1d of diastolic pressure in patients at risk of SCD (HRJSDds_P- HRJSDds_ 0.0014 [0.0010; 0.0018] 0.0053 [0.0035; 0.0071] *** LD1-E2 LU1). Similar patterns were also observed in the cardio-diastolic HRJSDcd_ 0.0004 [0.00028; 0.00053] 0.00014 [0.00010; 0.00019] *** coupling (HRJSDcd_LD1-P), suggesting that the deterioration E0-LA1 of autonomic regulation is more severe in patients at high risk HRJSDds_ 0.069 [0.063; 0.074] 0.093 [0.086; 0.10] *** of SCD. Earlier studies indicated that symmetric patterns in the LD1-LD1 HRJSD could be related to baroreflex-like response patterns. This suggests that this kind of behavior is also more pronounced in *ρ ≤ 0.01; **ρ ≤ 0.001; ***ρ ≤ 0.0000167; n.s. ρ ≥ 0.01; IDC, idiopathic cardiomyopathic patients; CON, control group; PPA, Poincaré plot analysis; HRJSD, high resolution patients at high risk (Baumert et al., 2005; Schulz et al., 2016). joint symbolic dynamic; NSTPDC, normalized short-time partial directed coherence; cd, The Poincaré plot analysis revealed that the patients from cardio-diastolic coupling; cs, cardio-systolic coupling; ds, diastolic-systolic coupling. the HR group have higher short-term systolic blood pressure deviation than the patients from the LR group and the CON Our findings revealed that the cardio-diastolic coupling subjects. This pathological behavior is also reflected in the (NSTPDCcd_NF) is bidirectional with the diastolic activity as PPAs_SD1/SD2 index, which is less balanced as the illness a driver in IDCLR patients. Additionally, the coupling strength progresses. The diastolic blood pressure behaved in a similar from cardiac activity over both the systolic and diastolic blood way: the short-term diastolic blood pressure deviation was pressure (NSTPDCcs_Area_c->s, NSTPDCcs_Area_c->d) gets significantly lower in HR patients and their PPAd_SD1/SD2 was stronger when the patients are at high risk, suggesting that the less balanced as well, indicating higher BPV in patients with cardiac activity is significantly dominating the blood pressure. critical conditions. In earlier studies (Tatasciore et al., 2013; This type of relationship has been observed before in congestive Ribeiro et al., 2017), higher short-term BPV was associated heart failure patients (Marinazzo et al., 2006). with several cardiac maladies such as left ventricular systolic Our results also revealed that there is a significant increase in dysfunction and atherosclerosis, amongst others. This BPV systolic pressure activity as a response to the alternant activity behavior was also present in sinoaortic denervated cats (Di Frontiers in Physiology | www.frontiersin.org 11 July 2019 | Volume 10 | Article 841 Rodriguez et al. Cardiovascular Coupling and Risk Stratification TABLE 8 | Mean value, 95% confidence interval and p-value of the HRV and BPV standard indices and the dual sequence method across all comparisons. Index IDCHR (N =14) IDCLR (N = 77) CON (N = 49) p-value p-value p-value p-value IDCHR vs. IDC vs. IDCHR vs. IDCLR vs. IDCLR CON CON CON bslope 5.36 [3.49; 7.22] 7.81 [6.8; 8.75] 10.29 [9.09; 11.49] n.s. ** ** * tslope 5.51 [3.75; 7.27] 7.94 [7.05; 8.83] 11.01 [9.78; 12.24] n.s. ** ** ** BBI_meanNN 828.42 [775.26; 881.57] 906.95 [879.27; 934.63] 883.79 [853.34; 914.23] n.s n.s n.s. n.s. BBI_sdNN 33.44 [25.06; 41.81] 36.98 [32.8; 41.09] 47.82 [43.45; 52.19] n.s. ** * ** BBI_rmssd 17.43 [12.00; 22.87] 20.53 [18.22; 22.84] 32.90 [28.25; 37.55] n.s. ** * ** BBI_pNN50 2.48 [1.45; 3.52] 3.44 [2.73; 4.14] 0.14 [0.09; 0.18] n.s. ** * ** BBI_HFn 0.65 [0.60; 0.70] 0.67 [0.63; 0.70] 0.62 [0.58; 0.67] n.s. n.s. n.s. n.s. BBI_LFn 0.34 [0.29; 0.39] 0.32 [0.29; 0.36] 0.37 [0.32; 0.41] n.s. n.s. n.s. n.s. BBI_LF/HF 2.48 [1.45; 3.52] 3.44 [2.73; 4.14] 2.55 [2.06; 3.04] n.s. n.s. n.s. n.s. SYS_meanNN 121.3 [107.6; 135] 112.6 [108.3; 117] 122.5 [116.7; 128.2] n.s. n.s. n.s. ** DIA_meanNN 61.5 [55.2; 67.9] 58.5 [55.8; 61.2] 60.6 [56.7; 64.5] n.s. n.s. n.s. n.s. DIA_VLF 0.48 [0.4;.0.56] 0.48 [0.45; 0.51] 0.39 [0.35; 0.44] n.s. ** n.s. ** *p ≤ 0.01; **p ≤ 0.001; n.s. p ≥ 0.001; IDCHR , idiopathic cardiomyopathic patients at high risk of sudden cardiac death; IDCLR , idiopathic cardiomyopathic patients at low risk of sudden cardiac death; CON, control group; bslope, bradycardic fluctuations; tslope, tachycardic fluctuations; BBI, beat-to-beat cardiac interval; SYS, systolic blood pressure; DIA, diastolic blood pressure. FIGURE 6 | SVM classification plots: (A) IDC vs. CON (Laplace kernel), (B) IDCHR vs. IDCLR (Laplace kernel), (C) IDCHR vs. CON (Gaussian kernel), and (D) IDCLR vs. CON (Laplace kernel). Rienzo et al., 1991). Additionally, baroreflex effectiveness has suggesting that the HRV in the CON group is higher than in IDC been studied in paraplegic subjects (Castiglioni et al., 2007), and it patients. Additionally, the BBI_rmssd index suggests that the their BPV, compared to control subjects, was found to be higher. HRV is lower in IDC patients compared with the CON subjects The aforementioned patterns were also present when the (Shaffer and Ginsberg, 2017). It is known that reduced HRV is a IDC patients were compared to the CON group. In general, predictor of an adverse prognosis in patients with cardiac disease the indices found in the PPA suggest that short-term BPV is (Task-Force-of-the-European-Society-of-Cardiology-the-North- higher in patients with pathological conditions (Mancia et al., American-Society-of-Pacing-Electrophysiology, 1996; Gang 1983; Mehlum et al., 2018). The PPA indices regarding the BBI and Malik, 2003). Several studies have related low HRV with revealed that short-term deviation and the short- and long-term heart failure (Dekker et al., 2000; Musialik-Lydka et al., 2003; deviation ratio of the heart rate is higher in the CON group, Sandercock and Brodie, 2006; Patel et al., 2017). Additionally, an Frontiers in Physiology | www.frontiersin.org 12 July 2019 | Volume 10 | Article 841 Rodriguez et al. Cardiovascular Coupling and Risk Stratification TABLE 9 | Accuracy (Acc), sensitivity (Sn), specificity (Sp), and area under the characterization methods are more suitable for describing the curve (AUC) obtained with the best SVM model for each comparison. non-linear behavior of HRV in pathological conditions. IDCHR vs. IDC vs. CON IDCHR vs. IDCHR vs. The HRJSD results suggest that patients at high risk adapt IDCLR CON CON less frequently to changes in blood pressure, reflected in the lower presence of decreasing patterns of BBI (E0) in response Indices PPAs_ HRJSDds- HRJSDds_ HRJSDds- to decreasing patterns in diastolic blood pressure (LD1). This SD1/SD2 E0d LD1-V E1d may be a result of the vagal response, causing less frequent HRJSDds_ NSTPDCcs_ NSTPDCcs_ NSTPDCcs_ LU1-P Area_c-s Area_c-s Area_c-s parasympathetic activity, leading to a less effective control of the C 7 5.5 1 0.3 blood pressure, and consequently the heart rate, in patients at Kernel Laplace Laplace Gaussian Laplace higher risk conditions. These results were consistent when CON subjects were σ 1 1.5 0.3 1 compared with IDC patients: decreasing patterns (E0) in Acc (%) 98.9 93.6 96.8 84.9 diastolic BPV were reflected in decreasing heart rates at higher Sn (%) 100 93.7 92.9 89.6 frequencies in CON subjects. The DIA_LF/HF was higher in Sp (%) 93.1 95.5 98.0 79.5 the patients when compared to the CON group. Higher levels AUC 0.96 0.94 0.95 0.85 of this index reflects efferent sympathetic activity (Mccraty and Shaffer, 2015). Additionally, a higher prevalence of unchanging patterns (E1) was found in the IDC patients compared with the CON group, indicating that changes in blood pressure are frequently not reflected in changes in heart rate in patients with pathological conditions. In addition, the HRJSDcd-E2d and HRJSDcd-LD1d indices showed that steady (E2) and low decreasing (LD1) diastolic blood pressure patterns, independent from all BBI patterns, are more frequent in CON subjects. These indices suggest a worsened circulatory homeostasis in IDC patients and support the idea of the influence of baroreflex activity in pathological conditions FIGURE 7 | Cascade model structure. The subject is evaluated in the IDC vs. (Kishi, 2018). CON model, and resulting IDC patients’ level of risk is evaluated by means of the HR vs. LR model. Earlier studies (Gavish et al., 2008; Schillaci and Pucci, 2010) have stated that the relationship between systolic and diastolic blood pressure should be coherent: if one of them increases, the other is expected to increase as well. The results increase complexity in BBI randomness and a lower fractal-like of the systolic-diastolic coupling revealed that patterns that are behavior have been associated with SDC (Huikuri et al., 2000). opposing in nature (sLU1-dP, sLD1-dV) are more frequent in the The DSM revealed that both the tslope and bslope are HR group. This suggests that the relationship between systolic significantly lower in IDC patients compared to CON subjects, and diastolic blood pressure loses linearity as the pathological and lower values of these sequences are associated with baroreflex condition worsens. dysfunction (Di Rienzo et al., 2001; La Rovere et al., 2008). In The coupling strength of the cardiovascular and diastolic- addition, a trend emerged in these indices when the patients were systolic couplings are stronger in pathological conditions. In compared by their level of risk: IDCHR patients showed lower addition, the symmetric patterns of the diastolic-systolic coupling tslopes and bslopes. activity were less frequent in patients. This may be caused by the The segmented approach of the Poincaré Plot analysis effect of the autonomic regulation mechanisms in pathological revealed that some patterns in the cardiovascular coupling conditions (Floras, 2009; Kishi, 2012). are more common in HR patients. There was a significantly To summarize, our results suggest that there is a gradual higher concentration of these patterns in column 5 and loss of HRV as SCD risk increases and, at the same time, BPV 8 in both cardio-diastolic and cardio-systolic couplings in increases alongside SCD risk. There is great controversy about patients at higher risk of SDC, indicating a lower variability the prognostic value of linear time and frequency domain HRV in their baroreflex activity compared with patients at low indices for risk stratification among this type of patient (Grimm risk. SPPA patterns (SPPAcd_Column_5, SPPAcd_Column_8) et al., 2003; Hohnloser et al., 2003; Minamihaba et al., 2003; occurred more frequently and were more concentrated in low Valencia et al., 2009; Voss et al., 2010, 2012a). The results of this risk patients, suggesting that these patients present a higher HRV work support the idea that the commonly used techniques for compared to patients at higher risk of SCD. The viability of HRV analysis of the time and frequency domain of HRV is not suitable as a reliable predictor of SCD in IDC patients was questioned for risk stratification. However, the combination of non-linear in an earlier study (Grimm et al., 2003), which did not support HRV analysis and linear as well as non-linear coupling analysis the hypothesis that HRV is a reliable predictor of SCD in IDC seems to be a promising tool for risk assessment in IDC patients patients. However, their results were based on a time domain (Voss et al., 2012b; Valencia et al., 2013; Fischer and Voss, analysis of HRV only, whereas the significant indices analyzed 2014). The processes involved in the circulatory homeostasis in our work come primarily from non-linear methods. These are by nature non-linear. Therefore, the differences between Frontiers in Physiology | www.frontiersin.org 13 July 2019 | Volume 10 | Article 841 Rodriguez et al. Cardiovascular Coupling and Risk Stratification the stages of this process can be more adequately revealed by of risk stratification in IDC patients. We have introduced a novel the quantification of the signal properties rather than by the cascade model that successfully classified subjects into different assessment of their magnitude. levels of SCD risk (CON, IDCLR , and IDCHR ). Further, this study We hypothesize that a dysfunction of the vagal activity, and allowed us to uncover, for the first time, some of the complex in general of the baroreflex mechanism as a whole, prevents interactions that take place within autonomic regulation, leading the body from maintaining circulatory homeostasis correctly. to a more accurate modeling, and interpretation of these This reduction in vagal activity and increase in the sympathetic processes in pathological conditions. Our findings revealed that influence exposes the cardiovascular system to frequent states there is a gradual decrease of HRV and an increase of BPV of stress that contribute to the worsening of the condition as the SCD risk increases. We conclude that patients at high over time. The gradual worsening of the heart rate and blood risk of SCD can no longer maintain circulatory homeostasis pressure variability’s in the different SDC risk stages considered consistently, leading to states of stress that worsens the here supports this assumption. This kind of impairment has condition over time. also been associated with other cardiac pathological conditions However, these results should be validated with a greater like ventricular fibrillation during myocardial ischemia (La number of patients, especially in the high-risk group. Therefore, Rovere et al., 2001). Abnormal sympathetic neural firing has the results presented in this work are more of a hypothesis- been associated with SCD and the genesis of ventricular generating nature than confirmatory. tachyarrhythmias (Esler et al., 1995). A similar behavior has been observed in elder mice (Freeling and Li, 2015), who ETHICS STATEMENT presented a reduced baroreflex bradycardic response compared to younger mice. This study are part of German Autonomic Regulation Trial This study has some limitations that are important to (ART) study, oriented to evaluate risk predictors of SCD and consider. The average age in the CON group was lower than in to improve the risk stratification in IDC patients. Signals the IDC group, the influence of age in the study of HRV has been from 220 IDC patients were recorded in two hospitals: the widely studied in earlier work (Voss et al., 2015). However, this Friedrich-Schiller University Hospital in Jena (44 patients) limitation does not affect our results for HR vs. LR comparisons. and the Franz-Volhard Clinic in Berlin (176 patients). Each A higher number of indices were analyzed than the number of study was performed following the recommendations of the patients in the database. Consequently, in order to minimize Declaration of Helsinki and patient consent was granted. problems due to possible overfitting, we stablished different levels The study was approved by the local Ethics Committee of statistical significance to reduce the dataset dimensionality, (No. 0986-11/02). including the Bonferroni-Holm correction criterion. The characterization of linear and non-linear coupling could AUTHOR CONTRIBUTIONS be also analyzed based on linear and non-linear Granger causality in time- and frequency domains (Schulz et al., 2013a). The BG and AV: conceptualization. JR, SS, BG, and AV: methodology, objective of this study was to evaluate the general behavior of validation, writing—original draft preparation, writing—review the underlying coupling, through the average of all features and editing, and final version approval. BG, JR, and SS: software (windows) over time was performed. Nevertheless, it is uncertain and data analysis. AV: data acquisition. if another time-invariant time domain measure based on Granger causality would be more appropriate and would have more FUNDING discriminative power (Porta and Faes, 2016). It would be interesting for an ongoing study. This work was partly supported by grants from the German Another limitation of this work is related to comorbidities Federal Ministry for Economic Affairs and Energy (KF and confounding factors influencing the autonomic regulation 2447309KJ4, ZF4485201SB7), by the Secretariat of Universities system. Therefore, these exclusion criterions make risk and Research of the Department of Economy and Knowledge stratification not yet applicable to every patient. of the Government of Catalonia (Consolidated research group GRC 2017 SGR 1770), by the CERCA Program/Generalitat CONCLUSION de Catalunya, and by the Spanish Ministry of Economy and Competitiveness (DPI2015-68820-R MINECO/FEDER), With this study we suggest that indices from coupling analysis and the Spanish Ministry of Science, Innovation and and non-linear HRV and BPV can contribute to the development Universities (RTI2018-098472-B-I00). REFERENCES Aro, A. L., Anttonen, O., Tikkanen, J. T., Junttila, M. J., Kerola, T., Rissanen, H. A., et al. (2011). Intraventricular conduction delay in a standard 12-lead Adochiei, F., Schulz, S., Edu, I., Costin, H., and Voss, A. (2013). A new normalised electrocardiogram as a predictor of mortality in the general population. short time PDC for dynamic coupling analyses. Biomed. Tech. 58(Suppl. 1). Circ. Arrhythm. Electrophysiol. 4, 704–710. doi: 10.1161/CIRCEP.111. doi: 10.1515/bmt-2013-4167 963561 Frontiers in Physiology | www.frontiersin.org 14 July 2019 | Volume 10 | Article 841 Rodriguez et al. Cardiovascular Coupling and Risk Stratification Aro, A. L., Anttonen, O., Tikkanen, J. T., Junttila, M. J., Kerola, human cardiac sympathetic nerves to stressors. Circulation 91, 351–358. T., Rissanen, H. A., et al. (2012a). Prevalence and prognostic doi: 10.1161/01.CIR.91.2.351 significance of T-wave inversions in right precordial leads of a 12-lead Faes, L., Erla, S., and Nollo, G. (2012). Measuring connectivity in linear electrocardiogram in the middle-aged subjects. Circulation 125, 2572–2577. multivariate processes: definitions, interpretation, and practical analysis. doi: 10.1161/CIRCULATIONAHA.112.098681 Comput. Math. Methods Med. 2012:140513. doi: 10.1155/2012/140513 Aro, A. L., Huikuri, H. V., Tikkanen, J. T., Junttila, M. J., Rissanen, H. Faes, L., Porta, A., Nollo, G., and Javorka, M. (2017). Information decomposition A., Reunanen, A., et al. (2012b). QRS-T angle as a predictor of sudden in multivariate systems: definitions, implementation and application to cardiac death in a middle-aged general population. Europace 14, 872–876. cardiovascular Networks 19, 1–28. doi: 10.3390/e19010005 doi: 10.1093/europace/eur393 Fischer, C., and Voss, A. (2014). Three-dimensional segmented poincaré Baccala, L. A., and Sameshima, K. (2001). Partial directed coherence: a plot analyses SPPA3 investigates cardiovascular and cardiorespiratory new concept in neural structure determination. Biol. Cybern. 84, 463–474. couplings in hypertensive pregnancy disorders. Front. Bioeng. Biotechnol. 2:51. doi: 10.1007/PL00007990 doi: 10.3389/fbioe.2014.00051 Baumert, M., Baier, V., Truebner, S., Schirdewan, A., and Voss, A. (2005). Floras, J. S. (2009). Sympathetic nervous system activation in human heart failure: Short- and long-term joint symbolic dynamics of heart rate and blood clinical implications of an updated model. J. Am. Coll. Cardiol. 54, 375–385. pressure in dilated cardiomyopathy. IEEE Trans. Biomed. Eng. 52, 2112–2115. doi: 10.1016/j.jacc.2009.03.061 doi: 10.1109/TBME.2005.857636 Freeling, J. L., and Li, Y. (2015). Age-related attenuation of parasympathetic Baumert, M., Walther, T., Hopfe, J., Stepan, H., Faber, R., and Voss, A. (2002). Joint control of the heart in mice. Int. J. Physiol. Pathophysiol. Pharmacol. 7, 126–135. symbolic dynamic analysis of beat-to-beat interactions of heart rate and systolic Fujita, H., Acharya, R., Sudarshan, V., Ghista, D., Sree, V., Lim Wei Jie, L., blood pressure in normal pregnancy. Med. Biol. Eng. Comput. 40, 241–245. et al. (2016). Sudden cardiac death (SDC) prediction based on nonlinear doi: 10.1007/BF02348131 heart rate variability features and SCD index. Appl. Soft Comput. 43, 210–519. Ben-Hur, A., Ong, C. S., Sonnenburg, S., Scholkopf, B., and Ratsch, G. (2008). doi: 10.1016/j.asoc.2016.02.049 Support vector machines and kernels for computational biology. PLoS Comput Gang, Y., and Malik, M. (2003). Heart rate variability analysis in general medicine. Biol. 4:e1000173. doi: 10.1371/journal.pcbi.1000173 Indian Pacing Electrophysiol. J. 3, 34–40. Bertini, M., Schalij, M. J., Bax, J. J., and Delgado, V. (2012). Emerging Gavish, B., Ben-Dov, I. Z., and Bursztyn, M. (2008). Linear relationship between role of multimodality imaging to evaluate patients at risk for systolic and diastolic blood pressure monitored over 24 h: assessment and sudden cardiac death. Circ. Cardiovasc. Imaging 5, 525–535. correlates. J. Hypertens. 26, 199–209. doi: 10.1097/HJH.0b013e3282f25b5a doi: 10.1161/CIRCIMAGING.110.961532 Giraldo, B. F., Rodriguez, J., Caminal, P., Bayes-Genis, A., and Voss, A. (2015). Bertinieri, G., Di Rienzo, M., Cavallazzi, A., Ferrari, A. U., Pedotti, A., and Cardiorespiratory and cardiovascular interactions in cardiomyopathy patients Mancia, G. (1988). Evaluation of baroreceptor reflex by blood pressure using joint symbolic dynamic analysis. Conf. Proc. IEEE Eng. Med. Biol. Soc. monitoring in unanesthetized cats. Am. J. Physiol. 254, H377–H383. 2015, 306–309. doi: 10.1109/EMBC.2015.7318361 doi: 10.1152/ajpheart.1988.254.2.H377 Grimm, W., Herzum, I., Muller, H. H., and Christ, M. (2003). Value of Bristow, M. R., Saxon, L. A., Boehmer, J., Krueger, S., Kass, D. A., De Marco, T., heart rate variability to predict ventricular arrhythmias in recipients of et al. (2004). Cardiac-resynchronization therapy with or without an implantable prophylactic defibrillators with idiopathic dilated cardiomyopathy. Pacing Clin. defibrillator in advanced chronic heart failure. N. Engl. J. Med. 350, 2140–2150. Electrophysiol. 26, 411–415. doi: 10.1046/j.1460-9592.2003.00060.x doi: 10.1056/NEJMoa032423 Hohnloser, S. H., Klingenheben, T., Bloomfield, D., Dabbous, O., and Cohen, Castiglioni, P., Di Rienzo, M., Veicsteinas, A., Parati, G., and Merati, G. (2007). R. J. (2003). Usefulness of microvolt T-wave alternans for prediction of Mechanisms of blood pressure and heart rate variability: an insight from low- ventricular tachyarrhythmic events in patients with dilated cardiomyopathy: level paraplegia. Am. J. Physiol. Regul. Integr. Comp. Physiol. 292, R1502–1509. results from a prospective observational study. J. Am. Coll. Cardiol. 41, doi: 10.1152/ajpregu.00273.2006 2220–2224. doi: 10.1016/S0735-1097(03)00467-4 Chugh, S. S. (2017). Sudden cardiac death in 2017: spotlight on prediction and Huikuri, H. V., Makikallio, T. H., Peng, C. K., Goldberger, A. L., Hintze, U., and prevention. Int. J. Cardiol. 237, 2–5. doi: 10.1016/j.ijcard.2017.03.086 Moller, M. (2000). Fractal correlation properties of R-R interval dynamics and Cohen, M. A., and Taylor, J. A. (2002). Short-term cardiovascular oscillations mortality in patients with depressed left ventricular function after an acute in man: measuring and modelling the physiologies. J. Physiol. 542, 669–683. myocardial infarction. Circulation 101, 47–53. doi: 10.1161/01.CIR.101.1.47 doi: 10.1113/jphysiol.2002.017483 Javorka, M., Krohova, J., Czippelova, B., Turianikova, Z., Lazarova, Z., Javorka, Cortes, C., and Vapnik, V. (1995). Support vector networks. Mach. Learn. 20, K., et al. (2017). Basic cardiovascular variability signals: mutual directed 273–297. doi: 10.1007/BF00994018 interactions explored in the information domain. Physiol. Meas. 38, 877–894. Dekker, J. M., Crow, R. S., Folsom, A. R., Hannan, P. J., Liao, D., Swenne, doi: 10.1088/1361-6579/aa5b77 C. A., et al. (2000). Low heart rate variability in a 2-minute rhythm strip Kishi, T. (2012). Heart failure as an autonomic nervous system dysfunction. J. predicts risk of coronary heart disease and mortality from several causes: the Cardiol. 59, 117–122. doi: 10.1016/j.jjcc.2011.12.006 ARIC Study. Atherosclerosis risk in communities. Circulation 102, 1239–1244. Kishi, T. (2018). Disruption of central antioxidant property of nuclear factor doi: 10.1161/01.CIR.102.11.1239 erythroid 2-related factor 2 worsens circulatory homeostasis with baroreflex Di Rienzo, M., Parati, G., Castiglioni, P., Omboni, S., Ferrari, A. U., Ramirez, dysfunction in heart failure. Int. J. Mol. Sci. 19:646. doi: 10.3390/ijms19030646 A. J., et al. (1991). Role of sinoaortic afferents in modulating BP and pulse- Kober, L., Thune, J. J., Nielsen, J. C., Haarbo, J., Videbaek, L., Korup, E., et al. interval spectral characteristics in unanesthetized cats. Am. J. Physiol. 261, (2016). Defibrillator implantation in patients with nonischemic systolic heart H1811–H1818. doi: 10.1152/ajpheart.1991.261.6.H1811 failure. N. Engl. J. Med. 375, 1221–1230. doi: 10.1056/NEJMoa1608029 Di Rienzo, M., Parati, G., Castiglioni, P., Tordi, R., Mancia, G., and Pedotti, La Rovere, M. T., Bigger, J. T. Jr., Marcus, F. I., Mortara, A., and Schwartz, A. (2001). Baroreflex effectiveness index: an additional measure of P. J. (1998). Baroreflex sensitivity and heart-rate variability in prediction of baroreflex control of heart rate in daily life. Am. J. Physiol. Regul. total cardiac mortality after myocardial infarction. ATRAMI (Autonomic Tone Integr. Comp. Physiol. 280, R744–R751. doi: 10.1152/ajpregu.2001.28 and Reflexes After Myocardial Infarction) Investigators. Lancet 351, 478–484. 0.3.R744 doi: 10.1016/S0140-6736(97)11144-8 Duray, G., Israel, C. W., and Hohnloser, S. H. (2006). Recent primary prevention La Rovere, M. T., Pinna, G. D., Hohnloser, S. H., Marcus, F. I., Mortara, A., implantable cardioverter defibrillator trials. Curr. Opin. Cardiol. 21, 15–19. Nohara, R., et al. (2001). Baroreflex sensitivity and heart rate variability in the doi: 10.1097/01.hco.0000198978.55637.b5 identification of patients at risk for life-threatening arrhythmias: implications Ebrahimzadeh, E., Pooyan, M., and Bijar, A. (2014). A novel approach to predict for clinical trials. Circulation 103, 2072–2077. doi: 10.1161/01.CIR.103.16.2072 sudden cardiac death (SCD) using nonlinear and time-frequency analyses from La Rovere, M. T., Pinna, G. D., Maestri, R., Mortara, A., Capomolla, S., Febo, HRV signals. PLoS ONE 9:e81896. doi: 10.1371/journal.pone.0081896 O., et al. (2003). Short-term heart rate variability strongly predicts sudden Esler, M. D., Thompson, J. M., Kaye, D. M., Turner, A. G., Jennings, G. L., cardiac death in chronic heart failure patients. Circulation 107, 565–570. Cox, H. S., et al. (1995). Effects of aging on the responsiveness of the doi: 10.1161/01.CIR.0000047275.25795.17 Frontiers in Physiology | www.frontiersin.org 15 July 2019 | Volume 10 | Article 841 Rodriguez et al. Cardiovascular Coupling and Risk Stratification La Rovere, M. T., Pinna, G. D., and Raczak, G. (2008). Baroreflex sensitivity: Roes, S. D., Borleffs, C. J., Van Der Geest, R. J., Westenberg, J. J., Marsan, measurement and clinical implications. Ann. Noninvasive Electrocardiol. 13, N. A., Kaandorp, T. A., et al. (2009). Infarct tissue heterogeneity 191–207. doi: 10.1111/j.1542-474X.2008.00219.x assessed with contrast-enhanced MRI predicts spontaneous ventricular Malberg, H., Wessel, N., Hasart, A., Osterziel, K. J., and Voss, A. (2002). Advanced arrhythmia in patients with ischemic cardiomyopathy and implantable analysis of spontaneous baroreflex sensitivity, blood pressure and heart rate cardioverter-defibrillator. Circ. Cardiovasc. Imaging 2, 183–190. variability in patients with dilated cardiomyopathy. Clin. Sci. 102, 465–473. doi: 10.1161/CIRCIMAGING.108.826529 doi: 10.1042/cs1020465 Rundle, J. B., Giguere, A., Turcotte, D. L., Crutchfield, J. P., and Donnellan, A. Malberg, H., Wessel, N., Schirdewan, A., Osterziel, K. J., and Voss, A. (1999). (2019). Global seismic nowcasting with shannon information entropy. Earth [Dual sequence method for analysis of spontaneous baroreceptor reflex Space Sci. 6, 191–197. doi: 10.1029/2018EA000464 sensitivity in patients with dilated cardiomyopathy]. Z. Kardiol. 88, 331–337. Sandercock, G. R., and Brodie, D. A. (2006). The role of heart rate variability in doi: 10.1007/s003920050294 prognosis for different modes of death in chronic heart failure. Pacing Clin. Mancia, G., Ferrari, A., Gregorini, L., Parati, G., Pomidossi, G., Bertinieri, G., Electrophysiol. 29, 892–904. doi: 10.1111/j.1540-8159.2006.00457.x et al. (1983). Blood pressure and heart rate variabilities in normotensive and Saour, B., J. Y., Philbin E, Tan, H., Nguyen, D., O’brien, J., Sidhu, S., et al. hypertensive human beings. Circ. Res. 53, 96–104. doi: 10.1161/01.RES.53.1.96 (2017). Following theraoeutic hypothermia in patients with sudden cardiac Marinazzo, D., Pellicoro, M., and Stramaglia, S. (2006). Nonlinear parametric death due to ventricular tachycardia or fibrillation. Int. J. Clin. Med. 8, 293–305. model for Granger causality of time series. Phys. Rev. E Stat. Nonlin. Soft Matter doi: 10.4236/ijcm.2017.85028 Phys. 73:066216. doi: 10.1103/PhysRevE.73.066216 Schillaci, G., and Pucci, G. (2010). The dynamic relationship between systolic and Mccraty, R., and Shaffer, F. (2015). Heart rate variability: new perspectives on diastolic blood pressure: yet another marker of vascular aging? Hypertens. Res. physiological mechanisms, assessment of self-regulatory capacity, and health 33, 659–661. doi: 10.1038/hr.2010.95 risk. Glob. Adv. Health Med. 4, 46–61. doi: 10.7453/gahmj.2014.073 Schmidt, A., Azevedo, C. F., Cheng, A., Gupta, S. N., Bluemke, D. A., Mehlum, M. H., Liestol, K., Kjeldsen, S. E., Julius, S., Hua, T. A., Rothwell, P. M., Foo, T. K., et al. (2007). Infarct tissue heterogeneity by magnetic et al. (2018). Blood pressure variability and risk of cardiovascular events and resonance imaging identifies enhanced cardiac arrhythmia susceptibility death in patients with hypertension and different baseline risks. Eur. Heart J. in patients with left ventricular dysfunction. Circulation 115, 2006–2014. 39, 2243–2251. doi: 10.1093/eurheartj/ehx760 doi: 10.1161/CIRCULATIONAHA.106.653568 Milde, T., Schwab, K., Walther, M., Eiselt, M., Schelenz, C., Voss, A., Schmidt, G., Malik, M., Barthel, P., Schneider, R., Ulm, K., Rolnitzky, L., et al. (2011). Time-variant partial directed coherence in analysis of the et al. (1999). Heart-rate turbulence after ventricular premature beats as a cardiovascular system. A methodological study. Physiol. Meas. 32, 1787–1805. predictor of mortality after acute myocardial infarction. Lancet 353, 1390–1396. doi: 10.1088/0967-3334/32/11/S06 doi: 10.1016/S0140-6736(98)08428-1 Minamihaba, O., Yamaki, M., Tomoike, H., and Kubota, I. (2003). Severity Schneider, T., and Neumaier, A. (2001). Algorithm 808: ARfit—a in myocardial dysfunction contributed to long-term fluctuation of heart matlab package for the estimation of parameters and eigenmodes of rate, rather than short-term fluctuations. Ann. Noninvasive Electrocardiol. 8, multivariate autoregressive models. ACM Trans. Math. Softw. 27, 58–65. 132–138. doi: 10.1046/j.1542-474X.2003.08207.x doi: 10.1145/382043.382316 Mozaffarian, D., Benjamin, E. J., Go, A. S., Arnett, D. K., Blaha, M. J., Schulz, S., Adochiei, F. C., Edu, I. R., Schroeder, R., Costin, H., Bar, K. J., et al. Cushman, M., et al. (2015). Heart disease and stroke statistics−2015 update: (2013a). Cardiovascular and cardiorespiratory coupling analyses: a review. a report from the American Heart Association. Circulation 131, e29–322. Philos. Trans. A Math. Phys. Eng. Sci. 371:20120191. doi: 10.1098/rsta.2012.0191 doi: 10.1161/CIR.0000000000000152 Schulz, S., Bär, K. J., and Voss, A. (2015a). Analyses of heart rate, respiration and Musialik-Lydka, A., Sredniawa, B., and Pasyk, S. (2003). Heart rate variability in cardiorespiratory coupling in patients with schizophrenia. Entropy 17, 483–501. heart failure. Kardiol. Pol. 58, 10–16. doi: 10.3390/e17020483 Neumaier, A., and Schneider, T. (2001). Estimation of parameters and eigenmodes Schulz, S., Bolz, M., Bar, K. J., and Voss, A. (2016). Central- and autonomic of multivariate autoregressive models. ACM Trans. Math. Softw. 27, 27–57. nervous system coupling in schizophrenia. Philos. Trans. Math. Phys. Eng. Sci. doi: 10.1145/382043.382304 374:20150178. doi: 10.1098/rsta.2015.0178 Panikkath, R., Reinier, K., Uy-Evanado, A., Teodorescu, C., Hattenhauer, Schulz, S., Haueisen, J., Bar, K. J., and Andreas, V. (2015b). High-resolution joint J., Mariani, R., et al. (2011). Prolonged Tpeak-to-tend interval on the symbolic analysis to enhance classification of the cardiorespiratory system in resting ECG is associated with increased risk of sudden cardiac death. patients with schizophrenia and their relatives. Philos. Trans. Math. Phys. Eng. Circ. Arrhythm. Electrophysiol. 4, 441–447. doi: 10.1161/CIRCEP.11 Sci. 373:20140098. doi: 10.1098/rsta.2014.0098 0.960658 Schulz, S., Tupaika, N., Berger, S., Haueisen, J., Bär, K. J., and Voss, A. (2013b). Parati, G., Di Rienzo, M., and Mancia, G. (2000). How to measure baroreflex Cardiovascular coupling analysis with high-resolution joint symbolic dynamics sensitivity: from the cardiovascular laboratory to daily life. J. Hypertens. 18, in patients suffering from acute schizophrenia. Physiol. Meas. 34, 883–901. 7–19. doi: 10.1097/00004872-200018010-00003 doi: 10.1088/0967-3334/34/8/883 Patel, V. N., Pierce, B. R., Bodapati, R. K., Brown, D. L., Ives, D. G., and Stein, P. Schwartz, P. J., La Rovere, M. T., and Vanoli, E. (1992). Autonomic nervous system K. (2017). Association of holter-derived heart rate variability parameters with and sudden cardiac death. Experimental basis and clinical observations for the development of congestive heart failure in the cardiovascular health study. post-myocardial infarction risk stratification. Circulation 85, I77–91. JACC Heart Fail. 5, 423–431. doi: 10.1016/j.jchf.2016.12.015 Seeck, A., Baumert, M., Fischer, C., Khandoker, A., Faber, R., and Voss, A. Porta, A., and Faes, L. (2016). Wiener-granger causality in network physiology (2011). Advanced Poincaré plot analysis differentiates between hypertensive with applications to cardiovascular control and neuroscience. Proc. IEEE 104, pregnancy disorders. Physiol. Meas. 32, 1611–1622. doi: 10.1088/0967-3334/ 282–309. doi: 10.1109/JPROC.2015.2476824 32/10/009 Porta, A., Faes, L., Bari, V., Marchi, A., Bassani, T., Nollo, G., et al. (2014). Shaffer, F., and Ginsberg, J. P. (2017). An overview of heart rate variability metrics Effect of age on complexity and causality of the cardiovascular control: and norms. Front. Public Health 5:258. doi: 10.3389/fpubh.2017.00258 comparison between model-based and model-free approaches. PLoS ONE. Silvani, A., Grimaldi, D., Vandi, S., Barletta, G., Vetrugno, R., Provini, F., et al. 9:e89463. doi: 10.1371/journal.pone.0089463 (2008). Sleep-dependent changes in the coupling between heart period and Ribeiro, A. H., Lotufo, P. A., Fujita, A., Goulart, A. C., Chor, D., Mill, J. G., et al. blood pressure in human subjects. Am. J. Physiol. Regul. Integr. Comp. Physiol. (2017). Association between short-term systolic blood pressure variability and 294, R1686–R1692. doi: 10.1152/ajpregu.00756.2007 carotid intima-media thickness in ELSA-brasil baseline. Am. J. Hypertens. 30, Stitson, M., Gammerman, A., Vapnik, V., Vovk, V., Watkins, C., and Weston, 954–960. doi: 10.1093/ajh/hpx076 J. (1999). Support Vector Regression With ANOVA Decomposition Kernels. Rodriguez, J., Voss, A., Caminal, P., Bayes-Genis, A., and Giraldo, B. F. (2017). Cambridge, MA: MIT Press. Characterization and classification of patients with different levels of cardiac Task-Force-of-the-European-Society-of-Cardiology-the-North-American- death risk by using Poincaré plot analysis. Conf. Proc. IEEE Eng. Med. Biol. Soc. Society-of-Pacing-Electrophysiology (1996). Heart rate variability: standards 2017, 1332–1335. doi: 10.1109/EMBC.2017.8037078 of measurement, physiological interpretation and clinical use. Task force of the Frontiers in Physiology | www.frontiersin.org 16 July 2019 | Volume 10 | Article 841 Rodriguez et al. Cardiovascular Coupling and Risk Stratification european society of cardiology and the North American society of pacing and Voss, A., Schulz, S., Schroeder, R., Baumert, M., and Caminal, P. (2009). Methods electrophysiology. Circulation 93, 1043–1065. derived from nonlinear dynamics for analysing heart rate variability. Philos. Tatasciore, A., Zimarino, M., Tommasi, R., Renda, G., Schillaci, G., Parati, G., et al. Trans. Math. Phys. Eng. Sci. 367, 277–296. doi: 10.1098/rsta.2008.0232 (2013). Increased short-term blood pressure variability is associated with early Wessel, N., Voss, A., Malberg, H., Ziehmann, C., Voss, H., Schirdewan, A., left ventricular systolic dysfunction in newly diagnosed untreated hypertensive et al. (2000). Nonlinear analysis of complex phenomena in cardiological patients. J. Hypertens. 31, 1653–1661. doi: 10.1097/HJH.0b013e3283 data. Z. Herzschr. Elektrophys 11, 159–173. doi: 10.1007/s0039900 61e4a6 70035 Valencia, J. F., Vallverdu, M., Porta, A., Voss, A., Schroeder, R., Vazquez, R., et al. Wolf, M. M., Varigos, G. A., Hunt, D., and Sloman, J. G. (1978). Sinus arrhythmia (2013). Ischemic risk stratification by means of multivariate analysis of the in acute myocardial infarction. Med. J. Aust. 2, 52–53. heart rate variability. Physiol. Meas. 34, 325–338. doi: 10.1088/0967-3334/3 Wu, K. C., Gerstenblith, G., Guallar, E., Marine, J. E., Dalal, D., Cheng, A., et al. 4/3/325 (2012). Combined cardiac magnetic resonance imaging and C-reactive protein Valencia, J. F., Vallverdu, M., Schroeder, R., Voss, A., Vazquez, R., Bayes De Luna, levels identify a cohort at low risk for defibrillator firings and death. Circ. A., et al. (2009). Complexity of the short-term heart-rate variability. IEEE Eng. Cardiovasc. Imaging 5, 178–186. doi: 10.1161/CIRCIMAGING.111.968024 Med. Biol. Mag. 28, 72–78. doi: 10.1109/MEMB.2009.934621 Wu, L., Jiang, Z., Li, C., and Shu, M. (2014). Prediction of heart rate variability Voss, A., Fischer, C., Schroeder, R., Figulla, H. R., and Goernig, M. (2010). on cardiac sudden death in heart failure patients: a systematic review. Int. J. Segmented Poincaré plot analysis for risk stratification in patients with dilated Cardiol. 174, 857–860. doi: 10.1016/j.ijcard.2014.04.176 cardiomyopathy. Methods Inf. Med. 49, 511–515. doi: 10.3414/ME09-02-0050 Voss, A., Fischer, C., Schroeder, R., Figulla, H. R., and Goernig, M. (2012a). Conflict of Interest Statement: The authors declare that the research was Lagged segmented Poincaré plot analysis for risk stratification in patients conducted in the absence of any commercial or financial relationships that could with dilated cardiomyopathy. Med. Biol. Eng. Comput. 50, 727–736. be construed as a potential conflict of interest. doi: 10.1007/s11517-012-0925-5 Voss, A., Goernig, M., Schroeder, R., Truebner, S., Schirdewan, A., and Figulla, Copyright © 2019 Rodriguez, Schulz, Giraldo and Voss. This is an open-access article H. R. (2012b). Blood pressure variability as sign of autonomic imbalance in distributed under the terms of the Creative Commons Attribution License (CC BY). patients with idiopathic dilated cardiomyopathy. Pacing Clin. Electrophysiol. 35, The use, distribution or reproduction in other forums is permitted, provided the 471–479. doi: 10.1111/j.1540-8159.2011.03312.x original author(s) and the copyright owner(s) are credited and that the original Voss, A., Schroeder, R., Heitmann, A., Peters, A., and Perz, S. (2015). Short-term publication in this journal is cited, in accordance with accepted academic practice. heart rate variability–influence of gender and age in healthy subjects. PLoS ONE No use, distribution or reproduction is permitted which does not comply with these 10:e0118308. doi: 10.1371/journal.pone.0118308 terms. Frontiers in Physiology | www.frontiersin.org 17 July 2019 | Volume 10 | Article 841