Paper accepted for publication on IJRRAS: International Journal of Research and Review in Applied Sciences A New Method for Analysis of Heart Rate Variability, Asimmetry and BRS. Part I . Elio Conte School of Advanced International Studies on Applied Theoretical and non Linear Methodologies of Physics, BariItaly and Department of Neuroscience and Sense Organs, University Aldo Moro, Bari-Italy Abstract: In the present paper we give a new method for estimation and quantification of Heart Rate Variability (HRV) in the VLF,LF,HF bands using the basic concept of variability previously introduced. The method enables to quantify ANS modulation of R-R intervals. In the subsequent paper we will give detailed exposition of the performed and confirming experiments. Key words : HRV, Analysis of ANS in HRV, New method of HRV analysis . Introduction The variation of the time period between consecutive heart beats, is thought to reflect the heart's adaptability to the changing physiological conditions. It is dependent predominantly on the extrinsic regulation of the heart rate [1]. Heart rate variability (HRV) analysis represents presently a powerful non-invasive method for analyzing the function of the autonomic nervous system. It is useful to understand the interplay between the sympathetic and parasympathetic autonomic nervous system, which serves to speed up and slow down the heart rate, respectively [1]. Assessment of HRV provides quantitative information about the modulation of heart rate (HR) by sympathetic nervous system (SNS) and parasympathetic nervous system (PNS). Interactions of SNS and PNS using HRV have been well studied and their importance established for a number of cardiac diseases including myocardial infarction [2], patients with congestive heart failure [3], patients at risk of sudden cardiac death [4,5] and patients with hypertension [6]. HRV analysis is currently used also to asses the dynamics of the HRV [7] describing the sympathetic and parasympathetic modulation of heart rate in other various clinical settings like diabetes [8], chronic heart failure [9] as previously said , chronic renal failure [10] and sleep apnea syndrome [11]. Heart rate variability analysis has valuable relevance also in psychophysiology and , in particular, in studies of physical, emotional, and mental function, as well as in psychological and psychiatric disorders[12]. In the present paper we take in consideration short terms (5-6 minutes) recording of beat-to beat intervals. There are two main standard approaches to the analysis of HRV: time-domain and frequency-domain analysis. Time-domain indices are derived from simple statistical calculations based on inter-beat intervals (RR intervals)). These indices are sensitive to transients and trends in the sample of heartbeats, and as such provide estimates of overall and beat-to-beat variability [13]. Let us outline here the most important indexes. The first index of Heart Rate Variability, widely used in medical research, is the Standard Deviation of the N-to-N intervals. The SDNN is the standard deviation of these intervals, a measure of their variability, expressed in milliseconds (ms). The RMSSD in msec is the square root of the mean of the sum of the squares of differences between adjacent NN intervals. Another measure of heart rate variability is the difference between the highest heart rate and the lowest heart rate within each cardiac cycle, measured in beats per minute. This index is called (HR Max – HR Min.) . Finally, a third index of variability, more reliable in the short term analysis , is called pNNX (usually pNN50). This index measures percent of difference between adjacent NN intervals and NN intervals that are greater than X(50) ms . Poincaré plot is still one of the popular time domain HRV analysis techniques. It is a visual presentation of time series signal to recognize the hidden patterns. The Poincaré plot of RR signal is constructed by plotting consecutive points of RR interval time series (i.e., lag-1 plot). It is a representation of the signal on phase space [14]. The most used technique available to quantify the Poincaré plot is fitting an ellipse to the shape of the Poincaré plot and measure the dispersion along the major and minor axis of the ellipse. By this way we realize a quantitative technique using two parameters , the short-term variability SD1 and the long-term variability SD2 to quantify the information from the plot. The technique was first proposed by Tulppo et. al. [15] . They defined two standard descriptors of the plot , SD1 and SD2 , for quantification of the Poincaré plot geometry. Later, the description of SD1 and SD2 in terms of linear statistics, given by Brennan et. al. [16], showed that the standard descriptors guide the visual inspection of the distribution. In case of HRV, it reveals a useful visual pattern of the RR interval data by representing both short and long term variations of the signal . . Frequencydomain analysis is currently the most employed technique in RR analysis. It is based on estimation of the power spectral density (PSD) of the RR time series by using FFT or DFT . The technique highlights the issue of the underlying rhythms of the mechanisms modulating heart rate and usually identifies three major spectral bands (high frequency (HF: 0.15-0.4 Hz), low frequency (LF: 0.04-0.15 Hz) and very low frequency (VLF: below 0.04 Hz)) in the adult spectrum . These measurements can be derived from short-term (i.e 5 to 6 minutes) ECG recordings. Long-term ECG (Holther) recordings (i.e. 24 hours) enables us to identify a further ULF (ultra low frequency band). The set of the previously described time and frequency techniques is currently used as a non-invasive marker of the activity of the autonomic nervous system . The necessary guidelines for comparing different studies of HRV have been established by the Task force of ESC and NPSE [17]. Some limitations may be identified in the current use of such methodologies. The RR signal is highly non linear, non periodic and non stationary [18] . The use of the non linear methodologies is required to place side by side standard linear methodologies. By using non linear methodologies , first of all a proper reconstruction of RR time series in phase space is required following the standard methodologies of the non linear analysis [18] . This procedure enables us to identify the embedding dimension first establishing a proper time delay by using Average Mutual Information and embedding dimension by using the FFN (False Nearest Neigh-boors) technique in order to account correct underlying temporal dynamics, the basic number of variables characterizing the dynamics of the system and , finally, the essential features of the reconstructed attractor. The primary limitation of the standard descriptors used for quantifying Poincaré plot is the lack of embedding temporal information. The standard descriptors, SD1 and SD2, represent the distribution of signal in 2D space and evaluate only information of width and length. In this manner also R-R signals having totally different underlying temporal dynamics may give similar Poincaré plots . Phase space reconstruction of the given RR time series usually enables to identify the Correlation Dimension ( )2D and the Lyapunov spectrum. In chaos theory, the correlation dimension is a measure of the dimensionality of the space occupied by the set of the given RR time series points. It may be referred to as a type of fractal dimension. The Lyapunov exponent or Lyapunov characteristic exponent of a dynamical system , as represented by the given RR time series, is a quantity that characterizes the rate of separation of infinitesimally close trajectories in the attractor . Quantitatively, two trajectories in phase space with initial separation zδ are considered to diverge at a rate given by zttz δλδ )exp()( = where λ is the Lyapunov exponent. The rate of separation can be different for different orientations of the initial separation vector. Thus, there is a spectrum of Lyapunov exponents that is equal in number to the dimensionality of the phase space. It is common to refer to the largest one as the Maximal Lyapunov exponent (MLE), because it determines a notion of predictability for a dynamical system. A positive MLE is usually taken as an indication that the system is chaotic and estimates its sensitivity to initial conditions . Actually, an arbitrary initial separation vector will typically contain some component in the direction associated with the MLE, and because of the exponential growth rate, the effect of the other exponents will be obliterated over time. Such indexes may be currently studied in RRR time series analysis provided that the experimental data are sufficiently controlled by performing an analysis of non linearity by AMI and using surrogate data as null hypothesis. The necessary guidelines for using such methodologies have been established by the Task force of ESC and NPSE [17]. There is still another methodology that is of high valuable interest in HRV analysis since it is able to look into the inner structure of the given RR time series. It is denominated Recurrence Quantification Analysis (RQA) as it was introduced by Webber and Zbilut [19] following the previous plot technique elaborated by Eckmann [20]. It is a method of non linear analysis which quantifies the number and the duration of recurrences of a dynamical system represented by its state space trajectory. The limits when using the non linear methodologies in RR time series analysis and HRV is that clinicians require detailed estimation of autonomic nervous modulation ( in particular sympathetic and parasympathetic autonomic nervous activity ) and valuable information in the VLF,LF,HF bands as they operate presently by using FFT or DFT that of course may present some limitations.. In order to make up to such insufficiency in our previous papers published on Chaos , Solitons and Fractals , we introduced a new methodology [21] that is denominated the CZF method also available as software by the Nevrokard company specialized in HRV analysis software (aHRV) (http://www.nevrokard.eu/) . However , in these years we have verified that in the application of the CZF method some clinicians encounter often some difficulties in the elaboration of the results as well as in the final conceptualization of the method. In order to avoid such further difficulty ,the aim of the present paper is to submit a new methodology devoted to the analysis of the variability in R-R signals in the frequency domain , including in particular the analysis in the VLF,LF , and HF bands in the frequency domain. The difference respect to the standard CZF method is that the present one represents a simplified but more rigorous and well arranged version of the CZF method so that clinicians will be able to apply it without further difficulties. 1. Elaboration of the methodology. Our study was essentially based on the attempt to estimate variability of a given RR time series by using non linear approaches. The previous CZF method attempted to give an answer to estimation of a so fundamental parameter, as the variability , but using the following algorithm . In essence, the Hurst analysis [22] brings to light that some statistical properties of a given time series X(t) scale with an observed period of observation T and a time resolution μ The scaling is characterized by a well known exponent H that relates the long-term statistical dependence of the signal. In substance we may express the scaling behaviour of statistically significant properties of the signal. indicating by E the mean values and considering to analyze the q-order moments of the distribution of the increments q qq tXE tXtXE k )(( )()(( )( −+ = τ τ (1.1) The (1.1) represents the statistical time evolution of the given stochastic variable )(tX . For q = 2, we may consider to rewrite the (1.1) in the following oversimplified form [ ] 2)( 1 )()( )(2 1)( ∑ = + −= hn i ihi uXuXhn hγ (1.2) The (1.2) estimates the variogram that we take in consideration where )(hn is the number of pairs at lag distance h while )( iuX and )( hiuX + are the time sampled series values at times t and ht + , ,........, 21 uut = . , ,........3,2,1=h . . .. . .. The variogram is thus a statistical measure in the form [ ] γ2)()( =−+ uXhuXVar Let us consider )(tX to represent a time series ( Dt∈ that is a subset of R+). As general scheme, let us suppose that this signal )(tX is composed by the sum of a deterministic part, )(tμ , plus a stochastic part. According to the previous relations we have that that )(2)()( huXhuXVar γ=−+ Consequently we may use the fractal variance function, )(hγ , the generalized fractal dimension, dimD , by the following equation: dim)( DhCh =γ and the simple marginal density function for self-affine distributions, given in the following manner: 1)( −−= aahkahP with dimD = 1−a and k is the scale parameter. This was the essence of the CZF method that we introduced previously [ 21 ]. By using the function [ ] 2)( 1 )()( )(2 1)( ∑ = + −= hn i ihi uXuXhn hγ (1.3) we performed essentially a study of the variability of the RR , represented by the relation )()( ihi uXuX −+ (1.4) Transferring the analysis in the frequency domain , we were able to estimate the VLF,LF,HF bands as usually obtained by application of the standard FFT or DFT but this time holding fundamentally about the concept of variability induced on the RR time series of a subject by the modulation operated by the ANS and its basic PNS and SNS interacting components. By the experience essentially based on some years of applications by the clinicians we have ascertained that the CZF , also if with large exceptions, usually gives difficulties to the clinical physiologist when it is the moment to conceptualize and reassume the obtained results. Consequently we will present a more direct and an oversimplified version in the present paper properly adapted for HRV studies and clinical applications. Let us recapitulate the question ab inito. Physiological systems are intrinsically non linear and exhibit high complexity. In healthy subjects such systems function under condition of disequilibrium , that is to say far from thermodynamic equilibrium and in these conditions exhibit selforganizing capability . In pathological conditions and diseases or in relation to the age a lowering or a loss of self-organizing capability is realized. As correctly outlined by Costa (22,Costa 2005) the self-organizing capability is related to time asymmetry of the underlying processes. Time asymmetry is here intended as lowering or loss of symmetry, meaning that the distribution of signals is in principle unbalanced and this basic property if progressively lowering or lost in pathological condition or disease and in age. In order to exemplify we may outline a general framework Biological organisms act far from equilibrium and they behave as complex systems controlled by non linear dynamics. The Deterministic Chaos Theory developed mathematical methods which have demonstrated to be useful studying complexity in biological signals Recently , the analysis of the Entropy and of Time Irreversibility has received valuable reconsideration. We started to evidence and to study the relevance of such problem of time irreversibility , also giving appropriate indexes for evaluation using radioactive tracers , starting with the 1980 and 1985 (The degree of Time asymmetry: on a new parameter of biological matter for estimation by radioisotope techniques and use in nuclear medicine)[23] Entropy quantifies regularity in a system, so a more regular series will be more predictable and less complex and its entropy will be lower. Low entropy reflects a less adaptable system and this may be observed in aging and illness. The auto-organizing capacity in a live organism is related to the uni-directionality of energy flow throw its systems and to the irreversibility of their processes. Time irreversibility consists in the loss of soundness in the statistic properties of a signal when one reverses its reading along the time. Two asymmetric trajectories are shown and the asymmetry index is higher in the healthy systems than in pathologic or aged one. The asymmetry in heart rate variability becomes an evident and characteristic phenomenon that we may directly appreciate in the Poincaré plot of normal sinus rhythm. The plot evidences the unevenness in the distribution of the points above and below the line of identity. In RR studies as well as in some general cases of physiological or clinical interest ,we may estimate and quantify the instantaneous changes in the beat to beat heart rate. This is of course the basic concept of variability that we take in consideration in this paper. So far , very little work has been devoted to analyze such kind of variability , linked to the previously mentioned asymmetry, in heart rate analysis. Let us examine the manner in which the two concepts , variability and time asymmetry, correlate. Variability may be intended by the following formula RR(n) –RR(n-k) (1.5) Where RR(n) states for the RR value at the beat n and RR(n-k) states for the RR value at the beat (n-k). k is an integer that may assume values k=1,2 ,3,4 , 5......usually denominated by us as lags. If we consider k=1 , we are estimating the difference RR(n) –RR(n-1) between the n beat and its previous one. The same concept holds for k=2, k=3, and so on. Essentially variability in RR holds about the estimation of the difference between a given RR values and its previous one or two preceding one , three preceding one and so on. This concept has obviously its direct counterpart in physiology where in fact RR(n) represents the time at the beat n and the previous or the subsequent fluctuations respect to such value, characterize directly the modulating action exerted basically from the PNS and SNS of the ANS. Of course the variability RR(n) –RR(n-k) represents the core of the CZF method that, as previously said, was based by us on the study of the (1.3) that in fact includes the (1.5). In the present simplified form of the CZF , let us admit that an )(nRR beat –to – beat time series is given (n=1,2,3,4,5,.... are the beats ). We suggest to estimate the previous conceptually established instantaneous variability of the )(nRR time series in the following manner )(nVaR = )( ))1()()(( nRRofValueMean nRRnRRnRR −− (1.6) To be explicit . Given the starting )(nRR beat –to – beat time series (n=1,2,3,4,5,.... are the beats ), for each beat we estimate variability ( ),...)5(),4(),3(),2( VaRVaRVaRVaR taking the value of the )(nRR time series multiplied for the variability that is represented from the difference of the given value minus the previous one and finally such value is divided ( and thus normalized) by the calculated mean value of the given )(nRR time series. )(nVaR is given in milliseconds and it evaluates the instantaneous value of the )(nRR variability at each beat. In this manner we reach the objective to study the variability and the asymmetry of the signal. Let us sketch the basic foundations of time asymmetry. Let us start giving the concept of a time series to be time reversible. Given the time series RR(n) . Let us admit that it is stationary . For RR(n) recording about five minutes it is usually assumed that subjects maintain such condition. The series RR(n) is said to be reversible if RR(1) ,RR(2), RR(3),........, RR(n) = RR(n),........,RR(3),RR(2),RR(1). Otherwise it is time irreversible. A trivial visual inspection to an RR(n) of an healthy subject evidences that this is not the case . Original RR (msec.) 600 700 800 900 1000 1100 1200 0 50 100 150 200 250 300 Evidence of time asymmetry in the time inverted RR 600 700 800 900 1000 1100 1200 0 50 100 150 200 250 300 In addition to visual inspection we have to quantify the phenomenon and to relate it to variability. A theorem may help us [24]. Let RR(n) be a time reversible series. Then for every k=1,2,3,...., the distribution of Var(n.k) with Var(n,k) = RR(n)-RR(n-k) (1.7) is symmetric about the origin. In conclusion we have reached two important results. The first is that we have found the manner to inspect if the given RR(n) time series is time symmetric or asymmetric , and previously we have outlined the importance of such estimation in clinical studies. We have to reconstruct the distribution of Var (n.k) and quantify such distribution respect to the origin. The second important feature is that by a theorem we have evidenced that time asymmetry and variability result profoundly linked since the (1.7) perfectly coincides with the (1.5) previously introduced by us when conceptualizing variability and its adherence with our starting CZF method. This last consideration completes the exposition of our CZF method revisited. The aim is now to give indexes of quantification that may be useful in RR and HRV analysis. First let us complete our estimation about the time asymmetry. Still a theorem helps us. As previously said , the function Var(n,k) (k=1,2,3,...... variability estimated at different lags) has a symmetric distribution about the origin for a time reversible process. The theorem is as it follows [24]. For each k=1,2,3,...... (k, variability estimated at different lags) we have that 0)).(()( == knVarsenofValueMeanhk ωω (1.8) for every +∈Rω In other terms, given the series RR(N) we may reconstruct the series Var(n,k). Soon after we may estimate the values of the function )),(( knVarsen ω and the arising mean value. If such calculated mean value results equal to zero we have time reversibility. Otherwise the process is time irreversible. Obviously the different obtained mean values are sufficient to study cases of clinical interest and, in particular, to discriminate normal from pathological conditions or normal subjects in function of the age. As final stage of our elaboration we have also to indicate some other methodological features. Let us confine our elaboration to variability estimation of two subsequent RR intervals. An alternative of the previous relation is to study also )(nVar taking the modulus of the difference )1()( −− nRRnRR that still may give important information respect to the case in which the modulus is not taken in consideration. Let us look at the sign of [ ])1()( −− nRRnRR . We may have [ ])1()( −− nRRnRR >0 or [ ])1()( −− nRRnRR < 0 The first case corresponds to a deceleration of the beat that we symbolize by (d) The second case corresponds instead to an acceleration that we symbolize by (a). Therefore, we are in the condition to estimate that total variability in the case of deceleration and in the case of acceleration. We indicate them by )(dVar and )(aVar , respectively. We may also compute the final asymmetry by ∑ ∑ ∑ ∑ + − = )()( )()( aVaRdVaR aVaRdVaR Asym as well as the coefficient of variation in both ( d) and (a ) cases respectively. There is still another parameter of interest. We may calculate the number of points in which the cardiac beat acquires deceleration as well as the number of cases in which it acquires acceleration and the number of points in which variability is not evidenced. The complex of such information may be useful in order to frame normal respect to pathological subjects. Finally , a further investigation to study also )(nVar taking the modulus of the difference )1()( −− nRRnRR that still may give important information respect to the case in which the modulus is not taken in consideration. We have completed the exposition of the basic indexes of our method. Let us consider now that the basic aim is to evaluate variability in the frequency domain and, in particular, evaluating Total Variability and the bands VLF, LF, and HF that are currently used in clinical applications. In our case we may proceed to frequency domain reconstruction by using the method that was exposed in [21] for the CZF. If the traditional customer does not intend to abandon instead the FFT or the DFT analysis we have to proceed as it follows. After the reconstruction of the time series )(nVaR we may thus proceed to perform its Fourier analysis by the DFT and thus evaluating the variability in the three bands of interest VLF,LF,HF , the total variability, the normalized unities and still more estimated ratios as LF/H F , LF/(LF+HF), HF(LF+HF) , collecting in this manner indexes that may be directly compared with those arising in standard procedure of the RR analysis. The basic difference is that using standard )(nRR time series , as we currently do, the estimated PSD results related to the variance of the )(nRR . In the case of the present method the variability is the key characterized by the given function )(nVaR = )( ))1()()(( nRRofValueMean nRRnRRnRR −− . And we estimate it in the VLF ,LF , and HF bands in order to evaluate the modulation as due to the autonomic nervous system (ANS). Let us introduce two figures to indicate the obtained results in the case of an health subject , 22 years and a 72 old subject considering the examined case without absolute value. It is clearly seen that the differences in Heart Rate Variability result strongly discriminated in relation to the age of the subject. In the upper left we have the starting tachogram resampled at 2 Hz. Down to the left we have )(dVar (deceleration) and )(aVar (acceleration) respect to the line of zero variability. In the upper right we have the DFT for the given tachogram. Instead , in the lower right we have the DFT for the variability , )(nVar . One observes the different values in amplitudes but the absolute coincidence in the peak location in Hz. Of course, one may perform time asymmetry analysis using one time the original tachogram and the other time the time inverted tachogram. This last statement concludes our exposition on the present, revisited and oversimplified version of our CZF method. We have now to give some comments on the manner to estimate baroreflex sensitivity (BRS). It is well known that the baroreflex loop is an important cardiovascular control mechanism for short-term blood pressure (BP) regulation. Based on afferent information of arterial baroreceptors reacting on changes in BP, central cardiovascular control is exerted on different peripheral effector systems as in particular on heart rate, cardiac output, peripheral resistance in order to keep BP between narrow limits. Baroreflex sensitivity (BRS) is a sensitive integrated measure of both sympathetic and parasympathetic autonomic regulation in which changes in heart rate due to variation in BP are reflected. Different techniques, based usually on spectral analysis [26] or on the so called sequence method [27], or correlation method [28] have been introduced to quantify baroreflex gain (29 ,Parlow et al. 1995). Traditionally, BRS is assessed pharmacologically, using the heart rate response to vasoactive drugs. Pharmacological and non invasive BRS measurements have been found to correlate significantly (30,Parlow et al. 1995, Watkins et al. 1996). The literature is boundless on this subject. However, no definitive agreement has been reached on which of the employed methods should be preferred (31,Lipman et al. 2003, Parati et al. 2004). On the other hand, BRS measurements represent an important prognostic tool to detect early subclinical autonomic dysfunction. In fact, reduced values of BRS may result largely from vagal withdrawal determining an high component of risk. Therefore it represents a valuable predictor of future cardiovascular morbidity and mortality in a variety of disease states that are associated with autonomic failure, (32,Gerritsen et al. 2001, La Rovere et al. 1998, Ditto 1990). An interesting connection has been shown between diminished BRS and psychopathology (33,Virtanen et al. 2003). We would suggest here the use of a method , based on the CZF [34] that accounts for variability in RR and SBP time series.. . Let us sketch briefly the essence of the approach. First we reconstruct variability and its representation in frequency domain as previously described in detail. The approach is as it follows : one time we will use our method on )(nRR time series and thus the same approach will be used for reconstructed SBP time series. In this condition, we will introduce two new BRS indexes of variability, BRSLF and BRSHF , for the two bands (LF) and (HF) respectively , calculated in the following manner: SBPforbandLFtheinyVariabilitTotal RRforbandLFtheinyVariabilitTotalBRSLF = in mmHg msec and SBPforbandHFtheinyVariabilitTotal RRforbandHFtheinyVariabilitTotalBRSHF = in mmHg msec A more accurate investigation requires obviously to use the mean blood pressure time series given by ((SBP − DBP) / 3 + DBP) instead of SBP. By using our method we have also the possibility to estimate the coupling strength of the variability existing simultaneously between RR and SBP . We may perform calculation by using only one lag as well as extending the investigation to include different lags. The procedure is as it follows. Given two time series, one for RR and the other for SBP , (we indicate them here by x (t) i and y (t) i ), the evaluation of the coupling strength and its sign will be estimate in the following manner [34, Conte ,Zbilut 2010] : Coupling Coefficient )()( )()( hh hhC SBPRR SBPRR SBPRR γγ γ − − = with [ ][ ][ ])()()()( 2 1)( tyhtytxhtxEhSBPRR −+−+=−γ The part II of the present paper will be devoted to the examination of the experimental results obtained by application of the present method in normal as well as in pathological subjects.