Multi- and monofractal indices of short-term heart rate variability R. Fischer I M. A k a y 2 P. Castiglioni 3 M. Di Rienzo 3 1Department of Biomedical Engineering, Rutgers University, Piscataway, USA 2Thayer School of Engineering, Dartmouth College, Hanover, USA 3LaRC, Centro di Bioingegneria, Fondazione Don C. Gnocchi, Milan, Italy Abstract--Indices of heart rate variability (HRV) based on fractal signal models have recently been shown to possess value as predictors of mortality in specific patient populations. To develop more powerful clinical indices of HRV based on a fractal signal model, the study investigated two HRV indices based on a monofractal signal m o d e l called fractional Brownian motion and an index based on a multifractal signal m o d e l called multifractional Brownian motion. The performance of the indices was compared with an HRV index in c o m m o n clinical use. To compare the indices, 18 normal subjects were subjected to postural changes, and the indices were compared on their ability to respond to the resulting autonomic events in HRV recordings. The magnitude of the response to postural change (normalised by the measurement variability) was assessed by analysis of variance and multiple comparison testing. Four HRV indices were investigated for this study: the standard deviation of all normal R-R intervals; an HRV index c o m m o n l y used in the clinic; detrended fluctuation analysis, an HRV index found to be the most powerful predictor of mortality in a study of patients with depressed left ventricular function; an HRV index developed using the m a x i m u m likelihood estimation (MLE) technique for a monofractal signal model; and an HRV index developed for the analysis of multifractional Brownian motion signals. The HRV index based on the MLE technique was found to respond most strongly to the induced postural changes (95% CI). The magnitude of its response (normalised by the measurement variability) was at least 25% greater than any of the other indices tested. Keywords--Heart rate variability, Detrended fluctuation analysis, Fractional Brownian motion, M a x i m u m likelihood estimation, Multifractal \ Med. Biol. Eng. Comput., 2003, 41,543-549 1 Introduction VARIABILITY IN the time interval between normal heartbeats provides information about the modulation of autonomic outflow to the heart. Mathematical metrics designed to characterise this variability have been shown to provide clinically significant information regarding the status of autonomic cardiac control and to be useful as a prognostic tool in a number of clinical settings. The analysis of heart rate variability (HRV) begins with a calculation of the time intervals between successive R-R intervals in a standard ECG record. The variability of the intervals is then characterised, typically with some time- or frequency-domain metric, in the time domain, for example, the variability can be calculated with the standard deviation of normal-to-normal (SDNN) interbeat time intervals, in the frequency domain, the variability is often calculated as the total spectral power within a predefined frequency band. Correspondence should be addressed to Professor Metin Akay; emaih metin.akay@dartmouth.edu Paper received 17 October 2002 and in final form 8 April 2003 MBEC online number: 20033800 © IFMBE: 2003 Medical & Biological Engineering & Computing 2003, Vol. 41 J in population studies, decreased HRV has been shown to possess predictive value for mortality in both the general population(TSUJI et al., 1994; 1996) and subjects with cardiovascular pathology. A relationship between decreased HRV and mortality in post-myocardial infarction (post-Mi) patients has been established in multiple studies (COPIE e t al., 1996; KLEIGERet al., 1987; VAISHNAVet al., 1994). Moreover, in the findings of the Multicenter Post-Infarction Project, HRV (using the SDNN index) was found to have the strongest univariate correlation with mortality of any Holter ECG variable (KLEIGERe t al., 1987). Despite the strong predictive value of decreased HRV, it may not have sufficient sensitivity and specificity to be used alone (TASK FORCEOF EUROPEAN SOCIETYOF CARDIOLOGY AND NORTH AMERICAN SOCIETYOF PACING ~; ELECTROPHYSIOLOGY,1996). in clinical work, an index of HRV is often coupled with other risk factors to improve its clinical utility. However, it may be possible to improve the sensitivity and specificity of HRV assessment by developing more powerful indices of HRV. This is a goal of the research presented here. 1.1 Characteristics of the HRV signal The value of the HRV measurement has inspired research efforts to associate aspects of the measurement with underlying 543 autonomic influences. Studies of the spectral content of HRV identified two major harmonic components: one at frequencies greater than 0.15 Hz, originally thought to be associated with the modulation of parasympathetic influences on HRV, and a component at or below 0.1 Hz, orignally thought to be associated with modulation of sympathetic influences on HRV (MALLIANI et al., 1991; SAUL et al., 1987). Broadband spectral power not associated with these two components was typically disregarded or believed to be measurement noise. This interpretation of the spectral content of HRV was reassessed however, in light of the discovery that HRV possessed a strong 1/f component over a wide frequency range KOBAYASHIand MUSHA, 1982). The 1/f component is present over many orders of magnitude, from time periods of 10 s to 10 h, and accounts for approximately 85% of the spectral power of resting human HRV (YAMAMOTO and HUGHSON, 1994). This finding has significant implications for researchers seeking to interpret features observed in HRV frequency spectra. First, the broad spectral aspects of the 1/fcomponent tend to bias the magnitude of peaks in HRV spectra, especially in the lowfrequency region. Secondly, the 1/f component implies a nonstationary process, meaning that simple HRV metrics such as SDNN will be influenced by the size of the data analysis window (YAMAMOTOand HUGHSON, 1994). Changes in the 1/f aspect of HRV apparently have physiological meaning as well. The LF coml~onent was found to be not strictly l/f, but more generally l/f", where fl is the spectral exponent, a parameter that determines the specific scaling behaviour in the frequency domain. Estimators of the spectral exponent were then derived, with the intention of providing useful information about the dynamics of the autonomic control of heart rate. Indeed, this is apparently the case; the spectral exponent of HRV has been shown to change in various physiological situations such as orthostatic challenge (BUTLER et al., 1993) ageing (LIPSITZe t al., 1990) and physical exercise (NAKAMURA et al., 1993). For example, experiments in which subjects were tested at increasing levels of head-up tilt displayed a progressive decrease in parasympathetic activity and a corresponding increase in spectral exponent (BUTLERe t al., 1993). 1.2 Fractal nature o f H R V The non-stationary character of HRV and its scaling aspect indicate that ffactal signal models may be appropriate for describing and analysing HRV. Typically, in a fractal signal analysis, we are interested in the calculation of certain scaling exponents that fully characterise the scaling properties. There are two general types of fractal signal model, monofractal and multifractal. 1.2.1 Monofractal signals: A monofractal signal can be characterised with a single scaling exponent. One of the most widely used monofractal signal models is the fractional Brownian motion (FBM) model. The FBM is a statistical signal model that includes many of the observed characteristics of HRV, including both the non-stationary aspect and scaling behaviour in the time and frequency domains. To characterise the variability of an FBM process, only a single parameter is required, the global Hurst exponent H. if we consider an HRV signal as the manifestation of an FBM process, we can then use the parameter H as the HRV index. indices of HRV based on the FBM model have been shown to carry prognostic information about mortality not extractable from traditional methods of HRV analysis in clinical studies of at-risk populations (MiKIKALLIO et al., 1997; 1998; 1999). in particular, an HRV index that uses a technique called detrended fluctuation analysis (DFA) to calculate H was 544 shown to be the best univariate predictor of mortality in a study of 159 patients with depressed left ventricular function following acute myocardial infarction (MAKIKALLIDe t al., 1999). 1.2.2 Multifractal signals: Multifractal signals have a highly complex structure that requires multiple scaling exponents for a complete description. A multifractal signal can be decomposed into many subsets, each characterised by different local Hurst exponents h. The goal of a multifractal analysis then is to calculate both the local exponents of the subsets and the fractal dimension associated with the organisation of the subsets. This information is described with a multifractal spectrum D(h), where D(h) is the fractal dimension of the subset of the signal having local Hurst dimension h. A common way to perform the analysis is by the wavelet transform. In a multifractal formalism developed by ARNEODOet al. (1991), signal singularities and the associated scaling exponents are recovered from the scaling of the wavelet transform along the so-called modulus maxima line. The modulus maxima line connects local maxima in the wavelet transform across the multiple scales of the wavelet decomposition. The scaling of the partition function Z(s, q) can be estimated from the modulus maxima values. The partition function is defined as the sum of the qth powers of the set of local modulus maxima at scale s. For small scales, the partition function has an approximate relationship with the scaling function z(q) on the moments q Z(s, q) ~_ s ~(q) (1) By linear regression on a log-log plot, the scaling function is recovered and then used to derive the multifractal spectrum through a Legendre transform D(h) = qh - z(q) (2) For healthy subjects, it was shown (IVANOVet al., 1999) that a multiffactal spectrum calculated in this manner on long (~6 h) records of HRV includes a broad range of local scaling exponents. However, for subjects with congestive heart failure, the range of scaling exponents is narrow, which is characteristic of a monofractal signal. These results imply that, under healthy conditions, HRV has a multiffactal nature and that a loss of multifractality may be associated with cardiac pathology. However, although this approach admits the presence of variability in the local Hurst exponent, it is statistical in nature and provides only a global estimate of scaling (STRUZlK,2000a). in some applications, however, the scaling properties can be non-stationary and can change in response to physiological events. For these applications, we would like information about local changes in scaling, and a single global estimate of scaling would not be useful. Fortunately, recent research has identified procedures that provide a stable local estimate of scaling and are also consistent with the global multifractal analysis approach (BENASSIe t al., 1998; STRUZlK, 2000a; b). With this in mind, our interest lies in analysis techniques that focus on local scaling information and permit investigation into the temporal variation of the scaling, and yet are consistent with a multifractal interpretation of HRV. Conventional multifractal analysis could not be used to investigate local scaling owing to the limited data available in local analysis data windows. However, other techniques are available that provide a stable estimate of the local scaling properties expected in a multifractal signal. In one approach, a generalisation of fractional Brownian motion called multifractional Brownian motion (MBM) has been developed (AYACHEand VEHEL, 2000; BENASSI et al., 1998). The major difference between FBM and MBM is that, for MBM, Medical & Biological Engineering & Computing 2003, Vol. 41 the scaling exponent can vary with time. However, an important attribute of the MBM model with regard to its analysis is that it has been shown to behave locally around a given point as a fractional Brownian motion (BENASSI et al., 1998; PELTIERand VEHEL, 1995). Thus FBM estimators can be used for the calculation of a local estimate of the Hurst exponent. 1.3 Research goals In the research presented here, our goal was to identify HRV indices that are consistent with a multifractal signal model for HRV and that are appropriate for the study of short-term autonomic influences on HRV. As stated previously, the use of FBM estimators for a local estimate of scaling is not inconsistent with a multiffactal interpretation of HRV. This is fortunate: powerful estimators of FBM have been developed that are appropriate for short time series (FISCHER and AKAY 1996; PILGRAM and KAPLAN, 1998). Fractal estimators can be classified into two categories, non-parametric and parametric. Non-parametric estimators typically depend on performing a linear regression fit of a scale-dependent quantity against scale in a logarithmic representation. For real-world fractals, the scaling relationship does not extend to arbitrarily large and small scales; consequently, estimators using this approach must select a relevant range of scales for analysis. Parametric estimators rely on a parametric representation of the autocorrelation function or spectral density. In this study, we utilised a non-parametric FBM estimator called detrended fluctuation analysis (DFA), an estimator found to be among the most powerful non-parametric methods in testing on simulated data (TAQQUe t at., 1995). We also investigated two parametric estimators. One estimator, developed using the MBM model, analyses the quadratic variations of a signal at a given point, and the other is a maximum likelihood estimator (MLE) for FBM. In the theory of statistical estimation, an MLE is often regarded as the best possible estimator, because it is asymptotically unbiased and efficient (HAYKIN, 1991). An important goal of our research was to determine whether the desirable properties of the MLE provided the basis for a useful HRV index in a clinical setting. The quantity produced by all these estimation techniques is h, the local value of the Hurst exponent, it is this value that is used as an index of heart rate variability. To test the effectiveness of these estimators as indices of HRV, we used data having known autonomic 'markers': clinical HRV recordings in which normal subjects underwent posmral shift at a specified time. it has been shown in previous studies that posmral tilt induces significant autonomic change that can be detected by HRV analysis (FEI et al., 1995). The autonomic change is thought to be the result of vagal withdrawal and increased sympathetic outflow (VYBIRAL et al., 1989). The HRV indices are then compared based on their response to this marker, relative to the inherent measurement variability. 2 Methods 2.1 H R V indices 2.1.1 Detrended fluctuation analysis estimator." The details of the DFA algorithm employed in this study have been described previously (PENG et al., 1995). The algorithm is a method for characterising the self-similar (or fractal) nature of a signal in the time domain. Briefly, in this algorithm, a time series is divided into segments of a given length. Each segment is 'detrended' by the subtraction of its linear trend. The root-mean-square fluctuation of each segment is then calculated and averaged over all segments. The average Medical & Biological Engineering & Computing 2003, Vol. 41 fluctuation is recalculated at a succession of segment lengths n and then plotted on a log-log scale against n. For a selfsimilar signal, this relationship will be linear, and the slope of the plot can be used to characterise the ffactal component and calculate a scaling exponent. The implementation of the DFA algorithm used here is the implementation developed by PENS et al. (1995), with one exception. The data submitted to the algorithm were differenced prior to analysis so that the scaling relationship was shifted by one. The calculated scaling exponent is then the Hurst exponent. 2.1.2 Maximum likelihood estimator." The theory of maximum likelihood estimation is derived from statistical estimation theory. In this approach, we considered each interbeat time interval to be the manifestation of a random process. The random process was represented by a probabilistic data model; this was the fractional Brownian motion model. As stated previously, the FBM model is dependent on a single parameter H. We used an MLE (described in Appendix 1) to estimate H given the input HRV data. 2.1.3 Multifractional Brownian motion estimator." The multifractal estimator employed in this study was developed based on a multifractional signal model called multifractional Brownian motion. The MBM model was developed independently by BENASSI et at. (1997) and PELTIER and LEVY-VEHEL (1995). The MBM model was developed to address a perceived shortcoming of the FBM model: the irregularity of an FBM is the same along its entire path. For an MBM, we no longer consider a global scaling exponent H for the entire path; instead, we calculate a local h that can vary with time. An estimator for MBM has been developed that calculates quadratic variations of the signal in the neighbourhood of a point t to produce an estimate of the scaling exponent, it is described in detail in Appendix 2. 2.1.4 Standard deviation o f normal interbeat intervals: The SDNN calculation was straightforwardly computed as the sample standard deviation of the RR intervals for each data window. 2.2 Data collection In this study, we used posmral changes to elicit changes in HRV that could be used as a basis for comparison of the HRV indices. Posmral change is a commonly used method to alter the balance between sympathetic and parasympathetic modulation of heart rate. The test subjects were asked to lie in a supine position, for approximately 10min in a quiet room with dim lighting, while electrocardiographic (ECG) signals were recorded. The subjects were then asked to sit upright as data were recorded for another 10 min interval. A posmral shift to an upright position is expected to elicit a change in autonomic outflow, such that parasympathetic modulation of HR is suppressed and sympathetic modulation is enhanced (VYBIRAL et al., 1989), increasing the value of h measured by the HRV indices. To assess the ability of the indices to respond to this posmral change, we used each estimator to calculate a series of HRV estimates before and after the autonomic event. The mean HRV estimates before and after the autonomic event could then be calculated to characterise the ability of the HRV index to respond to the underlying autonomic event. Although the magnitude of the index's response to the underlying autonomic event is important, it is also important to consider the response in the context of the inherent variability 545 supine upright = the same data set, the chance of erroneously declaring a result as statistically significant increases as we do more tests. The Bonferroni error protection method was used here to insure that the 'standard of proof' was made more stringent when we performed multiple comparison testing using a single data set (ZAR, 1996). 1.05E0.75~1.05 3 Results rr d- 0.75F 1.05~ 0.75~samplenumber Fig. 1 Typical HRV records for test subjects in upright and supine conditions of the measurement; in other words, to consider the signal-tonoise aspects of the HRV index. For a given HRV index, this was accomplished by calculating the following dimensionless ratio over the complete group of subjects: mean(HR V ,pright) SNR ---- - mean(HR V s,pi,e ) ~/var(HRV,pright) + var(HRVs,pi,e) (3) Data collection: Data were collected from 21 healthy, agematched subjects by ECG*, and finger arterial blood pressure was measured by a non-invasive device t. Each recording was digitised at 200 Hz and visually checked for artifacts. To calculate the RR intervals, a parabolic interpolation algorithm was used to refine the accuracy of the estimate of the time of occurrence of each R-wave. Abnormal QRS complexes were quite rare in the data recordings, occurring at a frequency of approximately 0.0004%. Records with a large proportion of artifacts were not considered; this resulted in three subjects being removed from this study. Preceding the HRV analysis, each of the 36 data records (18 subjects in two conditions: upright and supine) was windowed to a 512-point record length. This was the largest common window size that could be extracted from all data records. Typical R-R interval data for subjects in the upright and supine conditions are shown in Fig. 1. Statistical analysis: For each subject, data records before and after the postural shift were analysed with the four HRV indices. To compute mean and variance for the estimates produced by each index, a sequence of data analysis windows of 128-point length were used to calculate the HRV estimates. The analysis windows were overlapped by 50%, resulting in a total of 14 HRV estimates for each subject, seven before the postural shift and seven after. The estimates were then used to calculate mean and variance values for the estimates and to calculate SNR using (3). One-way analysis of variance was then performed on the signal-to-noise ratio (SNR) values for all subjects for the indices to determine if there was a significant difference in the SNR values produced by a particular HRV index. To determine which indices had higher SNR, comparisons of the DFA, MLE and MBM indices were performed using the SDNN index as control at a confidence interval of 95%. in this situation, however, it is important to adjust the error level to reflect the fact that, as multiple comparison tests are performed on hypotheses from *Cardioline Delta 1 Plus, Remco Italy tFinapres, Ohmeda 2300, Ohmeda Monitoring Englewood, CO, USA 546 The alteration in h induced by the postural change for all subjects is shown graphically in Fig. 2. Each line segment in the Figure represents the change in estimated HRV for a given subject due to the postural change. Although the true value ofh is unknown for these subjects, we can make some general observations. For the DFA and MLE indices, the value of h always increased for the supine to upright transition, with the exception of one subject. For the MBM index, however, the expected shift in h was reversed for this subject as well as two others. From the Figure, we can also compare the distribution of the h estimates between the three indices. Relative to the MLE index, the estimates produced by the DFA index are biased low. The estimates produced by the MBM index are biased high, beyond the expected range of 0.0-1.0. For statistical testing, the SNR was calculated for the 18 subjects for the indices and then analysed with the one-way analysis of variance. At ap-level of 0.014, the analysis found that the SNR of the indices was not equal. Multiple contrast testing (using the SDNN index as control) was then performed to determine which indices had significantly greater SNR. Using the Bonferroni error protection criterion at a confidence level of 95%, it was found that the MLE index responded more strongly to the postural changes than the SDNN index. At this confidence level, neither the DFA nor the MBM indices were found to differ significantly from the SDNN index. The statistical analysis is summarised in Table 1. 4 Discussion and conclusions The discovery that the HRV signal has a strong fractal component has inspired research into a variety of techniques for characterising the ffactal aspect. HRV indices based on these techniques have already achieved success in clinical applications. The success of these methods may lie in the fact that they examine more than just the magnitude of HRV variability (they characterise the correlation features of HRV at multiple time scales) and perhaps are a better way to characterise the status of the physiological control systems that are responsible for HRV. Studies of relatively long records of HRV ( ~ 6 h ) have demonstrated the multifractal nature of HRV. A complete //~ 0.8- 1.0- 0.60.40.20 Z//~ S supine upright supineupright supineupright Fig. 2 Changes in Hurst exponent H induced by postural shift from Systems, supine to upright condition, as estimated by (0) DFA, (•) MLE and (EB)MBM HRV indices Medical & Biological Engineering & Computing 2003, Vol. 41 Table 1 One-way ANOVA analysis o f HRVestimator response to postural change. Analysis finds that SNRs of various estimators" are not all equal (p = 0.014) and that only SNR for MLE estimator is significantly larger than SDNN estimator SNR by HRVmetric MLE DFA SDNN MBM n 18 18 18 18 Source of variation HRV estimator Within cells Total SSq 256 155 181 Contrast MLE against SDNN DFA against SDNN MBM against SDNN Difference 1.60 1.07 0.513 Mean 2.52 1.99 0.917 1.43 SD 2.35 1.53 0.716 0.848 SE 0.555 0.361 0.169 0.200 DF 3 68 71 MSq 8.67 2.28 F 3.80 P 0.0140 Bonferroni 95% CI 0.37-2.84 (significant) 0.16-2.31 0.72-1.75 description of HRV at this time scale requires a spectrum of scaling exponents. This spectrum reflects the fact that there is variability among the local scaling exponents and it describes the distribution of these exponents. In the research presented here, we have addressed the scaling behaviour of HRV on a much shorter time scale, in the order of a few minutes. Here, we are concerned with tracking changes in the local scaling exponent instead of describing the global distribution of scaling exponents. Our results show that an index based upon a monofractal signal model is best at detecting short-term alterations to autonomic balance. These results encourage further development of HRV indices based upon mono fractal signal models and highlight the power of such indices for the analysis of physiological events in HRV. Another outstanding question regards the accuracy o f a fractal model for HRV. Although effective indices of HRV can be derived from fractal models, it is by no means certain that a fractal model gives a complete description of the HRV signal. Certainly, HRV has other signal components, such as measurement noise, and harmonic components, such as respiratory sinus arrhythmia. Signal models that include these components may be the basis for more effective indices of HRV. Mathematically, to find the ML estimate, we need to express the probability density of the model as a function of the parameters (the likelihood function) and then attempt to find values of the parameters that maximise this function. By utilising the increments of FBM, we create a process that is zero-mean, normally distributed, strict-sense stationary and more amenable to analysis. The probability density function for the increments of FBM is apparently unimodal (LUNDAHL et al., 1986), and so, by taking the logarithm (a monotonic function) of the likelihood, we simplify our calculation without changing the location of the maximum. From (2), N 1 logp(x; H) = - ~ - l o g 2 n - ~loglRI - ~xTR 1X (4) with the elements of R available from the autocorrelation function as in (3). The two parameters to be estimated are H and o-2. if we decompose R as R = o-2R' (5) the likelihood function can be written as N Appendix 1 N 1 logp(x; H ) = - ~ - l o g 2 n - ~-log(o-2) - ~loglR '1 Maximum likelihood estimator f o r H We sought to derive an algorithm (the estimator) that would provide an estimate of H g i v e n the input HRV data. Further, we required that the algorithm have certain desirable statistical properties: estimates produced by an MLE can, under certain circumstances, be unbiased and have minimum variance. The derivation of the MLE starts with an expression for the autocorrelation function for the increments of FBM (LUNDAHL et al., 1986) r[k] v. = -~- Llk + 112H -- 2lkl 2H -[- Ik -- 112H] (1) With the autocorrelation known, the probability density function for the increments is straightforward (KAY, 1993) 1 { ~ p(x; H) = (2n)N/21RI1/2 exp -- } xrR 1x (2) 1 xrR,_l x (6) 2o-2 Equation (6) can be minimised over o-2 by taking the derivative with respect to o-2 and letting the derivative go to zero 82 xTR ~ 1x -- - - (7) N inserting this in (6) leaves us with the final expression of the likelihood that we will maximise over H N logp(x; H ) = - ~ - l o g 2 n - ~ - 1 2 loglR'l N 1 [xTR~ Ix~ og~ ~- -) N (8) where x is the dataset, R is the covariance matrix e i j = r [ l i - jl] (3) and r[k] is the autocorrelation function given in (1). Medical & Biological Engineering & Computing 2003, Vol. 41 To maximise this function over H, we need only search in the interval [0, 1]. For the implementation of the MLE algorithm used here, lower-upper (LU) decomposition of R was used to 547 find the inverse and determinant o f R, and Brent's algorithm was used to maximise the likelihood (PRESS et al., 1992). Appendix 2 Multifractional Brownian motion estimator The estimator calculates quadratic variations of the signal in the neighbourhood o f a point t to produce an estimate o f the scaling exponent. The signal X is observed at sampling points p / N , p = O, K, N. For any point o f time t during the observation o f the signal, a neighbourhood of t is defined as v~,N(t ) = { P - - t _< ~} (1) w h e r e p = 0 , . . . ,N, tE (0, 1), e > 0 , N > 0 . Local variations o f the signal are then calculated in the neighbourhood as the sum o f discrete second derivatives o f the signal Ve,N(t)= ~ (x(P+l,~_2x(P) + x(p~_))2 (2) This calculation is repeated again on a rescaled version o f the signal for which we consider only every other point, V~,N/Z (t). The local estimate o f the scaling exponent is then obtained as h(t) = 1 (log2 Ve'N/2(t) 1) -~ ~,N(t) ~- (3) Acknowledgem ents This work was supported by NIH grant (HL 65732). References ARNEODO, A., BACRY,E., and MuzY, J. (1991): 'Wavelets and multifractal formalism for singular signals: Application to turbulence data', Phys. Rev. Lett., 67, pp. 3515-3518 AYACHE, A., and VEHEL, J. (2000): 'The generalized multifractional Brownian motion', Stat. Infe~ Stochast. Process., 3, pp. 7-18 BENASSI, A., COHEN, S., and ISTAS, J. (1998): 'Identifying the multifractional function of a Gaussian process', Statist. Prob. Lett., 39, pp. 337-345 BENASSI, A., JAFFARD,S., and Roux, D. (1997): 'Gaussian processes and pseudodifferential elliptic operators', Revista Mathernatica iberoamericana, 13, pp. 19-89 BUTLER, G., YAMAMOTO, Y., XING, H., NORTHEY, D., and HUGHSON, R. (1993): 'Heart rate variability and fractal dimension during orthostatic challenges', J. Appl. Physiol., 75, pp. 2602-2612 CoPIg, X., HNATKOVA,K., STAUNTON,A., FEI, L., CAMM, A., and MALIK, M. (1996): 'Predictive power of increased heart rate versus depressed left ventriculax ejection fraction and heart rate variability for risk stratification after myocardial infarction. Results of a twoyear follow-up study', J Am. Coll. Cardiol., 27, pp. 270-276 FISCHER, R., and AKAY, M. (1996): 'A comparison of analytical methods for the study of fractional Brownian motion', Ann. Biomed. Eng., 24, pp. 537-543 HAYKIN, S. (1991): 'Adaptive filter theory, 2nd edn' (Prentice Hall, Englewood Cliffs, 1991) IVANOV, E, AMARAL, L., GOLDBERGER, A., HAVLIN, S., ROSENBLUM, M., STRUZIK, Z., et al. (1999): 'Multifractality in human heartbeat dynamics', Nature, 399, pp. 461-465 KAY, S. M. (1993): 'Fundamentals of statistical signal processing, vol. h estimation theory' (Prentice-Hall, Upper Saddle River, 1993) KLEIGER, R. E., MILLER, J. R, BIGGER, J. T., MOSS, A. J., and MULTICENTER POST-INFARCTION RESEARCH GROUP (1987): 548 'Decreased heart rate variability and its association with increased mortality after acute myocardial infarction', Am. J. Cardiol., 59, pp. 256-262 KOBAYASHI,M., and MUSHA, T. (1982): ' l / f fluctuation of heaxtbeat period', IEEE Trans. Biomed. Eng., 29, pp. 456-457 LIPSITZ, L. A., MIETUS, J., MOODY, G. B., and GOLDBERGER,A. L. (1990): 'Spectral characteristics of heart rate variability before and during postural tilt. Relations to aging and risk of syncope', Circulation, 81, pp. 1803-1810 FEI, L., ANDERSON, M. H., STATTERS, B. M., MALIK, M., and CAMM, A. J. (1995): 'Effects of passive tilt and submaximal exercise on spectral heart rate variability in ventriculax fibrillation patients without significant structural heart disease', Am. Heart J, 129, pp. 285-289 LUNDAHL, T., OHLEY, W. J., KAY, S. M., and SIFFERT, R. (1986): 'Fractional Brownian motion: A maximum likelihood estimator and its application to image texture', IEEE Trans. Med. Imag., 5, pp. 152-161 MALLIANI, A., PAGANI, E, LOMBARDI,E, and CERUTTI, S. (1991): 'Cardiovascular neural regulation explored in the frequency domain', Circulation, 84, pp. 482-492 MA.KIKALLIO, T. H., SEPPA.NEN, T., AIRAKSINEN, K. E. J., KOISTINEN, J., TULPPO,M. P., PENG, C. K., et al. (1997): 'Dynamic analysis of heart rate may predict subsequent ventricular tachycaxdia after myocaxdial infarction',Am. J. Cardiol., 80, pp. 779-783 M~.KIKALLIO,T. H., RISTIMA.E,T., AIRAKSINEN,K. E. J., PENG, C. K., GOLDBERGER, A. L., and HUIKURI, H. V (1998): 'Heart rate dynamics in patients with stable angina pectoris and utility of fractal and complexity measures', Am. J. Cardiol., 81, pp. 27-31 MA.KIKALLIO,T. H., HOIBER, S., KOBER, L., TORP-PEDERSEN, C., PENG, C. K., GOLDBERGER,A. L., et al. (1999): 'Fractal analysis of heart rate dynamics as a predictor of mortality in patients with depressed left ventriculax function after acute myocardial infarction', Am. J. Cardiol., 83, pp. 836-839 NAKAMURA,Y., YAMAMOTO,Y., and MURAOKA,I. (1993): 'Autonomic control of heart rate during physical exercise and fractal dimension of heaxt rate variability', J. Appl. Physiol., 74, pp. 875-881 PELTIER, R., and VEHEL, J. L. (1995): 'Multifractional Brownian motion: Definition and preliminary results'. Research Report 2645, INRIA National Institute for Research in Computer Science & Control, France PENG, C. K., HAVLIN,S., STANLEY,H. E., and GOLDBERGER,A. L. (1995): 'Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series', Chaos, 5, pp. 82-87 PILGRAM, B., and KAPLAN,D. T. (1998): 'A comparison of estimators for 1/f noise', Physica D, 114, pp. 108-122 PRESS, W. H., TEUKOLSKY, S. A., VETTERLING, W. T., and FLANNERY, B. E (1992): 'Numerical recipes in C, 2nd edn' (Cambridge University Press, New York, 1992) SAUL, J. R, ALBRECHT,E, BERGER, R. D., and COHEN, R. J. (1987): 'Analysis of long term heart rate variability: methods, 1/f scaling and implications', Cornput. Cardiol., 14, pp. 419-422 STRUZlK, Z. R. (2000a): 'Determining local singularity strengths and their spectra with the wavelet transform', Fractals, 8, pp. 163-179 STRUZlK, Z. R. (2000b): 'Revealing local variability properties of human heartbeat intervals with the local effective Holder exponent'. CWI Technical Report INS-R0015, National Research Institute for Mathematics & Computer Science, The Netherlands TAQQU, M. S., TEVEROVSKY, V, and WILLINGER, W. (1995): 'Estimators for long-range dependence: am emprical study', Fractals, 3, pp. 785-798 TASK FORCE OF EUROPEAN SOCIETYOF CARDIOLOGYAND NORTH AMERICAN SOCIETYOF PACING & ELECTROPHYSIOLOGY(1996): 'Heart rate variability. Standards of measurement, physiological interpretation, and clinical use', Circulation, 93, pp. 1043-1065 TsuJI, H., LAWRENCE, M. G., VENDITTI, E J., MANDERS, E. S., EVANS, J. C., FELDMAN,C. L., et al. (1994): 'Reduced heart rate variability and mortality risk in am elderly cohort. The Fraxningham Heart Study', Circulation, 90, pp. 878-883 TsuJI, H., LAWRENCE, M. G., VENDITTI, E J., MANDERS, E. S., EVANS, J. C., FELDMAN, C. L., et al. (1996): 'Impact of reduced heart rate variability on risk for cardiac events. The Fraxningham Heart Study', Circulation, 94, pp. 2850-2855 Medical & Biological Engineering & Computing 2003, Vol. 41 VAISHNAV, S., STEVENSON, R., MARCHANT, B., LAGI, K., RANJADAYALAN,K., and TIMMIS, A. (1994): 'Relation between heart rate variability early after acute myocardial infarction and long-term mortality', Am. J. Cardiol., 73, pp. 653-657 VYBIRAL, T., BRYG, R. J., MADDENS,M. E., and BODEN, W. E. (1989): 'Effect of passive tilt on sympathetic and parasympathetic components of heart rate variability in normal subjects', Am. J. CardioL, 63, pp. 1117-1120 YAMAMOTO,Y., and HUGHSON, R. L. (1994): 'On the fractal nature of heart rate variability in humans: effects of data length and fl-adrenergic blockade', Am. d Physiol., 35, pp. R40-49 ZAR, J. H. (1996). 'Biostatistical analysis, 3rd edn' (Prentice Hall, Upper Saddle River, 1996) Medical & Biological Engineering & Computing 2003, Vol. 41 Author's biography METIN AKAY is Associate Professor of Engineering, Psycology and Brain Sciences, and Computer Science at Dartmouth College, Hanover. He is author/coauthor of 14 books and is the founding and current editor-chief of the Biomedical Engineering Book Series published by the Wiley and IEEE Press. Dr. Akay's Neural Engineering and Informatics Lab is interested in investigating the respiratory somatosensory (RSS) responses of patients with obstructive sleep apnea syndrome (OSAS) and the effect of developmental abnormalities and maturation on the dynamics of respiration. He teaches courses in the areas of Biomedical Informatics, Signal Processing, Neural Networks and Computations. 549