Multiscale analysis of heart rate variability in non-stationary environments

Heart rate variability (HRV) is highly non-stationary, even if no perturbing influences can be identified during the recording of the data. The non-stationarity becomes more profound when HRV data are measured in intrinsically non-stationary environments, such as social stress. In general, HRV data measured in such situations are more difficult to analyze than those measured in constant environments. In this paper, we analyze HRV data measured during a social stress test using two multiscale approaches, the adaptive fractal analysis (AFA) and scale-dependent Lyapunov exponent (SDLE), for the purpose of uncovering differences in HRV between chronic fatigue syndrome (CFS) patients and their matched-controls. CFS is a debilitating, heterogeneous illness with no known biomarker. HRV has shown some promise recently as a non-invasive measure of subtle physiological disturbances and trauma that are otherwise difficult to assess. If the HRV in persons with CFS are significantly different from their healthy controls, then certain cardiac irregularities may constitute good candidate biomarkers for CFS. Our multiscale analyses show that there are notable differences in HRV between CFS and their matched controls before a social stress test, but these differences seem to diminish during the test. These analyses illustrate that the two employed multiscale approaches could be useful for the analysis of HRV measured in various environments, both stationary and non-stationary.

1. Introduction

Modern interest in heart rate variability (HRV) began when it was observed that it is more than an important and easily accessible indicator of cardiovascular function, but an important measure of autonomic nervous system (ANS) function and health in general. A salient feature is its spontaneous fluctuation, even if the environmental parameters are maintained constant and no perturbing influences can be identified. Obviously, variations in HRV will be more complicated if the data are measured in intrinsically non-stationary environments. Therefore new methods for better characterizing HRV measured in those situations are desirable.

Since the first observations that HRV can be a sensitive indicator of declining health that precedes changes in heart rate itself or other physiological measures of distress (Hon and Lee, 1963; Kleiger et al., 1987; King et al., 2009), a number of methods have been proposed to analyze HRV data. The most widely used methods assume a stationary process underlying the statistics, calculated from time and frequency domain analyses [see Malik (1996) and references therein], as well as those derived from chaos theory and random fractal theory (Kobayashi and Musha, 1982; Goldberger and West, 1987; Babyloyantz and Destexhe, 1988; Kaplan and Goldberger, 1991; Pincus and Viscarello, 1992; Bigger et al., 1996; Ho et al., 1997). In this paper, we illustrate the general use of two multiscale approaches that do not assume a stationary process, the adaptive fractal analysis (AFA) (Gao et al., 2011b; Kuznetsov et al., 2012; Riley et al., 2012), and the scale-dependent Lyapunov exponent (SDLE) (Gao et al., 2006a, 2007, 2012d), for the analysis of HRV in a study of chronic fatigue syndrome (CFS) patients with matched healthy controls.

Multiscale analysis of heart rate was first considered over a decade ago by groups exploring non-linear dynamics in physiology, for example (Costa et al., 2002, 2005). These studies, while quite interesting in their methods, explored datasets that were not subtle in their impact on cardiovascular physiology. Such sophisticated techniques were not needed to distinguish healthy heart rates from those in the diseases being considered (congestive heart failure and atrial fibrillation). The same dataset has been explored further by many others using different multiscale analysis techniques, for example (Hu et al., 2009a, 2010). These analyses have added to a growing interest in multiscale analysis of physiology in general (West, 2010), but studies using multiscale techniques to analyze more challenging datasets where the cardiac disturbances are more subtle (such as those found in CFS) are very rare. The techniques demonstrated here show some promise when the cardiac disturbances are subtle and the effect sizes are small.

CFS is a heterogeneous illness, with subsets of patients potentially having autonomic nervous system (ANS), immune system, and endocrine system involvement (Newton et al., 2011; Rahman et al., 2011). While no known biomarker for CFS has been identified to date, recent small-sample studies have suggested certain cardiac irregularities, such as low nocturnal HRV (Boneva et al., 2007; Burton et al., 2009; Rahman et al., 2011), small left ventricle and orthostatic intolerance (Miwa and Fujita, 2011), short QT interval (Naschitz et al., 2006), and low blood volume and diminished cardiac function (Hurwitz et al., 2009; Hollingsworth et al., 2011), are associated with CFS. It is not known if the cardiac irregularities are themselves a cause of the symptoms of CFS or if they are symptoms of the autonomic dysfunction and other systemic breakdowns which are more directly causative. Whatever the case, analysis of ECG data appears to show promise and, given the data accessibility and feasibility of collecting it, it behooves us to look beyond resting data to data collected in more stressful environments, and to analyze the data with methods that do not assume stationarity.

We shall employ two relatively new and fundamentally different multiscale approaches to characterize HRV of control and CFS patients. While we hope to shed some new light on the pathology of CFS through analysis of HRV, for the interest of this special issue, we will make our best efforts to illustrate how these methods are used to analyze our data, so that interested readers may readily adapt the methods to analyze their HRV data measured in other situations, both stationary and non-stationary.

2. Methods

2.1. Data

The HRV data analyzed here were derived from subjects who participated on the third day of the 3-day clinical study taking place at Emory University's Atlanta Clinical and Translational Science Institute (ACTSI, formerly known as General Clinical Research Center or GCRC). Subjects who participated in the ACTSI study were identified from the baseline (2004–2005) (Reeves et al., 2005) and follow-up waves of the Georgia longitudinal study on CFS and Chronic Unwellness. This study adhered to human experimentation guidelines of the U.S. Department of Health and Human Services, was approved by the Centers for Disease Control and Prevention (CDC) Human Subjects Review Board, and all participants gave written informed consent. The ACTSI study was a 1:2 case-control study design with matching on age within 5 years, sex, race/ethnicity, and BMI (≥30 kg/m², obese or not). The case-control status was determined by the classifications from the two waves (baseline and follow-up) of the Georgia longitudinal study. Subjects stayed overnight at the Emory clinic during their participation in the 3-day ACTSI study. The study consisted of two components: the first 2 days involved brain imaging in conjunction with studies of cognition and the ANS, and the third day involved a stress-induced challenge V the Trier Social Stress Test (TSST) V to test the hypothalamic pituitary adrenal (HPA) axis, ANS, and immune system. Between March 2008 and July 2009, there were 36 CFS cases and 48 non-fatigued healthy controls completed the 3-day study. The current analysis will focus on 23 CFS and 41 controls with the HRV data collected during the Day-3 TSST study. Heart rate data were collected using Biopac feeds into Somnologica software at 200 Hz on a Dell laptop computer. A technician attached sensors to subjects and calibrated equipment prior to the TSST. Subjects wore the electrodes and the recording device from 1:30 pm to 4:00 pm on Day 3. The ECG data were then stored and transferred to CDC computers. R–R interval pre-processing from the raw ECG was done with a LabVIEW script developed at the Brain-Body Center at University of Illinois, Chicago. This processing interpolated the R peaks from the 200 Hz data to improve their temporal localization.

Event times were annotated at the clinic in the Somnologica software itself, or, occasionally, in an Excel spreadsheet when the ECG recorder and/or Somnologica software was malfunctioning. Twelve 5-min intervals were chosen to represent each subject during the roughly 3 h of the TSST. These intervals correspond to three events (e.g., blood draws) before the test, four key events during the test (the receptionist's speech, a blood draw, the subject's TSST speech, and the math test), and five events after the test (blood draws), and allow correlations with gene expression to be made. All IBI data (inter-beat-interval, also called R–R interval, time between R peaks in msec) were manually inspected and corrected for artifacts prior to analyses (CardioEdit; Brain-Body Center, University of Illinois at Chicago, 2007). Missed and/or incorrect R peak detections were manually corrected using the ECG waveform when available. Otherwise, artifacts were corrected using integer arithmetic (i.e., dividing intervals when detections were missed and adding intervals when spuriously invalid detections occurred). During this process, event times were checked for discrepancies between the Somnologica annotations and the master Excel spreadsheet recording for all event times, and between those recorded times and other data (e.g., physiological events like rapid rise in heart rate during the TSST speech). Not all subjects had interpretable, cleanable raw ECG data, or well marked event times for all 12 of the data intervals. Some subjects suffered from arrhythmias which made the R peak picking software unreliable, and their data could not be used in this study. Furthermore, we excluded six subjects who were not consistently classified as CFS at some point during the multi-year study. Thus, in this study we used a reduced subset of all the subjects who underwent the TSST. Losses to the subject pool skewed somewhat the sex, age, race, and BMI matching of the study, but not significantly as can be seen in Table 1 below (p-values for age and BMI computed with a 2-sided t-test; those for sex and race using a 2-tailed Fisher's exact test).

TABLE 1

Table 1. Deomgraphics of subjects used in this study.

In this study, we will focus on the analysis of HRV data measured before and during the stress tests. To appreciate what our HRV data look like, we have shown in Figures 1A,B IBI data of a control and a CFS subject measured during the stress test. While the details of the two IBI data sets are different, we do observe some statistical similarity between them.

FIGURE 1

Figure 1. Examples of IBI data for a control (A) and CFS (B) subject. The red curves are the trend signals obtained by the adaptive filter, which will be explained below.

2.2. Adaptive Fractal Analysis (AFA)

In this paper we will characterize the fractal nature of the IBI time series with the Hurst parameter, H. Although many excellent methods for estimating H exist, including the famous detrended fluctuation analysis (Peng et al., 1994), care should be exercised by their interpretation, particularly if one is faced with relatively short time series that contain trends, non-stationarity, or signs of rhythmic activity (Hu et al., 2009b). One of the better approaches to tackle these difficulties is the AFA (Gao et al., 2011b, 2012a; Kuznetsov et al., 2012; Riley et al., 2012). In the following, we will attempt to explain the method clearly as we apply it to HRV data.

AFA is based on non-linear adaptive multiscale decomposition (Gao et al., 2010; Tung et al., 2011), which starts by partitioning a time series into segments of length w = 2n + 1, where neighboring segments overlap by n + 1 points, which ensures symmetry. Each segment is then fitted with the best polynomial of order M. Note that M = 0 and 1 correspond to piece-wise constant and linear fitting, respectively. We denote the fitted polynomials for the i-th and (i + 1)-th segments by y⁽ⁱ⁾ (l₁) and y^{(i + 1)} (l₂), respectively, where l₁, l₂ = 1, …, 2n + 1. We then define the fitting for the overlapped region as

$\begin{array}{cc}{y}^{(c)}(l)=w_1{y}^{(i)}(l+n)+w_2{y}^{(i+1)}(l),l=1,2,\dots ,n+1,& (1)\end{array}$

where $w_1=(1-\frac{l-1}{n})$ and $w_2=\frac{l-1}{n}$ can be written as (1 − d_j/n) for j = 1, 2, and where d_j denotes the distances between the point and the centers of y⁽ⁱ⁾ and y^{(i + 1)}, respectively. This means that the weights decrease linearly with the distance between the point and the center of the segment. Such a weighting ensures symmetry and effectively eliminates any jumps or discontinuities around the boundaries of neighboring segments. In fact, the scheme ensures that the fitting is continuous everywhere, is smooth at the non-boundary points, and has the right- and left-derivatives at the boundary. Moreover, since it can deal with an arbitrary trend without a priori knowledge, it can remove non-stationarity, including baseline drifts and motion artifacts (Gao et al., 2011b), and the procedure may also be used as either high-pass or low-pass filter with superior noise-removal properties than linear filters, wavelet shrinkage, or chaos-based noise reduction schemes (Tung et al., 2011). In Figures 1A,B we have plotted two trend signals in red for the IBI data shown there, using a window size of w = 101 and a polynomial order of 2.

Based on the described adaptive decomposition, a fractal analysis can be conducted as follows. Denote IBI data as u(i), i = 1, …, N, and the global smooth trend such as shown as red curves in Figures 1A,B as v(i), i = 1, …, N. AFA essentially is a scaling law relating the variance of the residual time series u(i)−v(i)and the window size w (Gao et al., 2011b)

$\begin{array}{cc}F(w)={\left[\frac{1}{N}{\displaystyle \sum _{i=1}^N{(u(i)-v(i))}^2}\right]}^{1/2}~w^H.& (2)\end{array}$

Note that in this formulation, IBI data are treated as random walk processes. This is consistent with the literature (Peng et al., 1993; Hu et al., 2010) and what is seen in Figures 1A,B. Also note that for truly fractal processes, the polynomial order does not matter. This is also largely true for IBI data. In this work, we always fix the polynomial order to be 1.

2.3. Scale-Dependent Lyapunov Exponent (SDLE) Analysis

SDLE is a multiscale complexity measure first introduced in 2006 (Gao et al., 2006a, 2007). it has been further developed theoretically (Gao et al., 2009, 2012b) and applied to characterize eeg (Gao et al., 2011a, 2012c), hrv (Hu et al., 2009a, 2010), financial time series (Gao et al., 2011c), and earth's geodynamo (Ryan and Sarson, 2008). recently, sdle is compared with a number of entropy measures (Gao et al., 2012c). it is found that sdle has superior scaling behaviors—it has well-defined scaling laws for all known major classes of time series, and embodies all the information approximate entropy and sample entropy may have. since we have done an in depth tutorial of sdle in Gao et al. (2012d), here we will only briefly describe sdle is such a way that the material presented here is self-contained.

SDLE is a concept defined in a high-dimensional phase space using the time delay embedding technique (Packard et al., 1980; Takens, 1981; Sauer et al., 1991). This technique is perhaps the most significant contribution of chaos theory to practical data analysis, since non-trivial dynamical systems usually involve many state variables, and therefore have to be described by a high-dimensional state (or phase) space. Consider a scalar time series x[n] = x(1), x(2), …, x(n). The embedding technique consists of creating vectors of the form:

$\begin{array}{cc}\begin{array}{l}{V}_i=[x(i),x(i+L),\dots ,x(i+(m-1)L)],i=1,\dots ,\\ N_p=n-(m-1)L,\end{array}& (3)\end{array}$

where n is the total length of the time series, N_p is the total number of constructed vectors, m is the embedding dimension and L the delay time. The embedding parameters need to be chosen according to certain optimization criteria. For details, we refer to chapter 13 of our book Gao et al. (2007).

After a proper phase space is re-constructed, we consider an ensemble of trajectories. We denote the initial separation between two nearby trajectories by ε₀, and their average separation at time t and t + Δt by ε_t and ε_{t + Δt}, respectively. We can then examine the relation between ε_t and ε_{t + Δt}, where Δt is small. When Δt → 0, we have,

$\begin{array}{cc}{\epsilon }_{t+\Delta t}={\epsilon }_t{e}^{\lambda ({\epsilon }_t)\Delta t},& (4)\end{array}$

where λ (ε_t) is the SDLE given by

$\begin{array}{cc}\lambda ({\epsilon }_t)=\frac{\ln {\epsilon }_{t+\Delta t}-\ln {\epsilon }_t}{\Delta t}.& (5)\end{array}$

Equivalently, we can express this as,

$\begin{array}{cc}\frac{d{\epsilon }_t}{dt}=\lambda ({\epsilon }_t){\epsilon }_t.& (6)\end{array}$

To compute SDLE, we check whether pairs of vectors (V_i, V_i) defined by Equation (3) satisfy the following Inequality,

$\begin{array}{cc}{\epsilon }_k\le \Vert V_i-V_j\Vert \le {\epsilon }_k+\Delta {\epsilon }_k,k=1,2,3,\dots ,& (7)\end{array}$

where ε_k and Δε_k are arbitrarily chosen small distances, and

$\begin{array}{cc}||V_i-V_j||=\sqrt{{\displaystyle \sum _{w=1}^m{\left(x_{i+(w-1)L}-x_{j+(w-1)L}\right)}^2}}& (8)\end{array}$

Geometrically, Inequality (Equation 7) defines a high-dimensional shell (which reduces to a ball with radius Δε_k when ε_k = 0; in a 2-D plane, a ball is a circle described by (x − a)² + (y − b)² = r², where (a, b) is the center of the circle, and r is the radius). We then monitor the evolution of all such vector pairs (V_i, V_j) within a shell and take the ensemble average over indices i, j. Since we are most interested in exponential or power-law functions, we assume that taking logarithm and averaging can be exchanged, then Equation (5) can be written as

$\begin{array}{cc}\lambda ({\epsilon }_t)=\frac{\langle \ln \Vert V_{i+t+\Delta t}-V_{j+t+\Delta t}\Vert -\ln \Vert V_{i+t}-V_{j+t}\Vert \rangle }{\Delta t}& (9)\end{array}$

where t and Δt are integers in units of the sampling time, the angle brackets denote the average over indices i, j within a shell, and

$\begin{array}{cc}{\epsilon }_t=\Vert V_{i+t}-V_{j+t}\Vert =\sqrt{{\displaystyle \sum _{w=1}^m{\left(x_{i+(w-1)L+t}-x_{j+(w-1)L+t}\right)}^2}}& (10)\end{array}$

To ease the following discussion, we now list two scaling laws for SDLE:

(1) For clean chaotic signal, λ(ε) fluctuates slightly around a constant. As is expected, this constant is the very largest positive Lyapunov exponent, λ₁,

$\begin{array}{cc}\lambda (\epsilon )={\lambda }_1.& (11)\end{array}$

(2) For noisy dynamics, on small scales,

$\begin{array}{cc}\lambda (\epsilon )~-\gamma \ln \epsilon ,& (12)\end{array}$

where γ is a coefficient controlling the speed of loss of information (i.e., defined as the measure of uncertainty involved in predicting the value of a random variable). This feature suggests that entropy generation is infinite when the scale ε approaches zero.

3. Results

3.1. AFA of HRV

Random fractal theory is perhaps the most famous model for HRV. A central theme of random fractal theory is the notion of long-range correlation characterized by the Hurst parameter H (Gao et al., 2006b, 2007): a time series is said to have anti-persistent, short-range or memoryless, or persistent long-range correlations if 0 < H < 1/2, H = 1/2, or 1/2 < H < 1, respectively. A classic result about HRV dynamics is that HRV data measured in quiet conditions possess anti-persistent correlations characterized by 0 < H < 1/2 (Peng et al., 1993).

Figure 2 shows two AFA curves for the data shown in Figures 1A,B. We observe that the scaling breaks around w = 2⁴. It turns out this is a generic feature among all the subjects. Note that the Hurst parameter H_s for the short-time scale is larger than 1/2, different from most of the literature that H for HRV is usually smaller than 1/2. The difference could be because our IBI data were measured during the TSST, while most published results (e.g., Peng et al., 1993; Hu et al., 2010) used HRV data measured in resting state. Another explanation could be that on short time scales HRV is more persistent, as is suggested by Poincare plots of heart rate, which typically show a strong positive correlation (Otzenberger et al., 1998; Smith and Reynolds, 2006; Smith et al., 2007). None-the-less, we observe that H_l is smaller than 1/2, indicating that on long time scales, IBI data of control and CFS subjects still have anti-persistent correlations, similar to HRV data measured in resting states. Finally, we note that the two AFA curves are very similar. This is because the original IBI data are similar, as we pointed out earlier.

FIGURE 2

Figure 2. AFA of the IBI data shown in Figures 1A,B. The scaling break around w = 2⁴ is generic among all the subjects. (A,B) The red circles are for computed H, the blue lines are a linear fit for H_s, and the green lines are linear fits for H_l.

The two very desirable properties of AFA are (1) it can readily deal with arbitrarily complicated trends, and (2) it works well even when the data set is short (Gao et al., 2012a). The latter property enables us to apply AFA to the six 5-min IBI data measured before and during the TSST (3 before and 3 during the TSST). When AFA is extended to these data sets, according to the two conditions, before and during the TSST, we find that the knee point around w≈ 2⁴ is still generically present. The second, long-time scale scaling regime, however, becomes shorter, since the IBI data are shorter. Nevertheless, H_s and H_l are still well-defined. The results are summarized in Figure 3, where (A1,A2) are for the data measured before the TSST, and (B1,B2) for the data measured during the stress tests. Simple statistical tests on the PDF's shown in Figure 3 are shown in Table 2. We observe that before the TSST, group differences exist for H_s and H_l between control and CFS subjects, although the area under the curve for the respective ROC plots indicates that each measure is only weakly discriminatory for individual subjects. However, neither H_s nor H_l appears to show any substantial differences between control and CFS during the stress tests (ROC curves were not computed during the TSST due to their low discriminatory power). These preliminary observations bear further investigation using studies with larger sample sizes.

FIGURE 3

Figure 3. Probability density functions (PDFs) for H_s and H_l of control and CFS subjects, where (A1, A2) are for data measured before the TSST, and (B1, B2) are for data measured during the stress tests.

TABLE 2

Table 2. Statistical tests on H_s and H_l before and during the TSST.

3.2. SDLE of HRV

In this section, we will focus on whether SDLE shows promise identifying differences between HRV of healthy control and CFS subjects during the TSST. For the purpose of the SDLE analysis, we ensured that 20–30 min of continuous IBI data during the TSST was obtained rather than only 5 min intervals. Such long data are needed for SDLE analysis. A drawback is that due to the small sample size, estimating the PDFs for metrics derived from SDLE analysis is ill-posed—this is more doable with AFA, since AFA works with 5-min IBI data, and we have several different data intervals before and during the TSST. Therefore with SDLE, we will use scatter plots.

Figure 4 shows two examples of error growth curves, ln ε(t) vs. t, and their corresponding SDLE curves, for a control and a CFS subject. These examples were chosen to illustrate some of the differences between subjects with regards to SDLE, and do not illustrate general differences between CFS cases and healthy controls. We observe that the general dynamics on short time scales is noisy dynamics characterized by a rapid increase in ln ε(t), and a scaling of Equation 12. Such behavior holds for all the subjects.

FIGURE 4

Figure 4. Error growth ln ε(t) vs. t and SDLE λ (ε) vs. ln ε curves for a control (A1, A2) and CFS (B1, B2) subject.

A number of metrics from the error growth and SDLE curves in Figure 4 can be derived and checked to determine if CFS and control groups can be separated in a statistically significant way. These metric are:

(1) Δε_max: It measures how far apart the error growth curves originating from different shells (or initial conditions) as shown in Figures 4A1,B1 are. It basically measured the degree of non-stationarity in the IBI data. The results are shown in Figure 5A.

(2) The coefficient γ₁ and γ₂ of Equation (12), for small and large ε scales, respectively. Here, small ε generally refers to the range of ε where the SDLE curve is almost straight. They are ln ε ≤ −2.6 and ln ε ≤ −3.2 for Figures 4A2,B2, respectively. The large ε scales for those plots correspond to ln ε ≥ −2.6 and ln ε ≥ −3.2. The results are shown in Figures 5B,C.

(3) ln ε_max: This corresponds to the largest value in the error growth curves, or the right most point in SDLE curves. Using an ensemble forecasting framework, we have proven that this quantify plays the role of the maximal amount of information (Gao et al., 2012c). The results are shown in Figure 5D.

(4) SDLE_max: the maximal SDLE value, shown in Figure 5E.

(5) SDLE_ε*: this is the SDLE at a fixed scale ε*. It plays a similar role as the commonly used approximate entropy and sample entropy (Gao et al., 2012c). The results are shown in Figure 5F.

FIGURE 5

Figure 5. (A–F) Summary of SDLE analysis, where blue circles and red stars are for control and CFS subjects, respectively (see text for more details about each metric).

In summary, for all the six metrics derived from SDLE analysis, there appears to be little difference between healthy control and CFS subjects. This appearance is easy to validate with simple statistical tests, as shown in Table 3 (p-values from a 2 tailed t-test).

TABLE 3

Table 3. Statistical tests on SDLE parameters during the TSST.

As can be seen in Table 3, none of the tests of the SDLE parameters to distinguish case from control groups are significant save one, and that would lose its significance after correction for multiple hypotheses using, e.g., a Bonferroni correction or a consideration of false discovery rate (FDR).

4. Discussion

CFS is a debilitating medical disorder withno known biomarker. Recent small-sample studies about possible cardiac irregularities (Naschitz et al., 2006; Boneva et al., 2007; Burton et al., 2009; Hurwitz et al., 2009; Hollingsworth et al., 2011; Miwa and Fujita, 2011; Newton et al., 2011; Rahman et al., 2011; Beaumont et al., 2012) have motivated us to carefully examine whether HRV in persons with CFS may be substantially different from that in healthy controls. Given that previous studies have already tried to quantify the difference between the HRV of healthy control and CFS subjects using standard HRV metrics that assume stationarity, we have chosen two completely different multiscale methods, AFA and SDLE, that do not assume stationarity to analyze HRV, and used them on data from healthy control and CFS subjects in a highly non-stationary environment. These two methods are fundamentally different, because AFA belongs to random fractal theory, while SDLE has its origin in deterministic chaos theory, although SDLE has been proven to also be able to characterize various types of random processes. Using AFA, the data suggests potential differences between the HRV of healthy control and CFS subjects prior to social stress tests, but the differences seem to be diminished during the test itself. The latter observation is further supported by SDLE analysis. Both observations may require further statistical validation and potentially larger sample sizes.

The differences observed in resting HRV between CFS and healthy controls (before the stress tests) is consonant with the observation that the complexity of HRV for CFS patients is reduced during sleep (Boneva et al., 2007). Furthermore, the differences in power spectral density of HRV seen by Boneva et al. (a greater loss of low frequency power in CFS vs. high frequency power) are consistent with a shifting of the Hurst parameter in CFS toward an anti-persistent regime seen in this study. While it might not be appropriate to conclude that there is no difference between the HRV of healthy control and CFS subjects during the TSST, the trend is clear—the difference between HRV of control and CFS subjects appears to be reduced during stress. This trend has also been suggested recently by other studies (Beaumont et al., 2012).

As we pointed out in the beginning, the purpose of our paper is twofold—to analyze HRV as a window into the potential ANS anomalies and pathophysiology of CFS, and to demonstrate the general use of AFA and SDLE for analyzing non-stationary HRV. The character of the stress test employed during this study naturally makes our HRV data more non-stationary than those measured in resting states. The results presented here suggest that AFA and SDLE could be very useful for the analysis of HRV, measured in both stationary and non-stationary environments.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors would like to greatly thank Dr. Stephen Porges of the University of Illinois at Chicago and the Research Triangle Institute for his help in reviewing and cleaning the data. We would also like to thank the CDC/Emory ACTSI Research Team for help with planning and collecting the data.

Disclaimer

The findings and conclusions in this report are those of the authors and do not necessarily represent the views of the CDC.

References

Babyloyantz, A., and Destexhe, A. (1988). Is the normal heart a periodic oscillator? Biol. Cybern. 58, 203–211. doi: 10.1007/BF00364139

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Beaumont, A., Burton, A., Lemon, J., Bennett, B., Lloyd, A., and Vollmer-Conna, U. (2012). Reduced cardiac vagal modulation impacts on cognitive performance in chronic fatigue syndrome. PLoS ONE 7:e49518. doi: 10.1371/journal.pone.0049518

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bigger, J., Steinman, R., Rolnitzky, L., Fleiss, J., Albrecht, P., and Cohen, R. (1996). Power law behavior of rr-interval variability in healthy middle-aged persons, patients with recent acute myocardial infarction, and patients with heart transplants. Circulation 93, 2142–2151. doi: 10.1161/01.CIR.93.12.2142

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Boneva, R., Decker, M., Maloney, E., Lin, J., Jones, J., Helgason, H., et al. (2007). Higher heart rate and reduced heart rate variability persist during sleep in chronic fatigue syndrome: a population-based study. Auton. Neurosci. Basic Clin. 137, 94–101. doi: 10.1016/j.autneu.2007.08.002

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Burton, C., Knoop, H., Popovic, N., Sharpe, M., and Bleijenberg, G. (2009). Reduced complexity of activity patterns in patients with chronic fatigue syndrome: a case control study. BioPsychoSocial Med. 3:7. doi: 10.1186/1751-0759-3-7

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Costa, M., Goldberger, A., and Peng, C. (2002). Multiscale entropy analysis of complex physiologic time series. Phys. Rev. Lett. 89:068102. doi: 10.1103/PhysRevLett.89.068102

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Costa, M., Goldberger, A., and Peng, C. (2005). Multiscale entropy analysis of biological signals. Phys. Rev. E 71:021906. doi: 10.1103/PhysRevE.71.021906

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gao, J., Cao, Y., Tung, W., and Hu, J. (2007). Multiscale Analysis of Complex Time Series—Integration of Chaos and Random Fractal Theory, and Beyond. Hoboken, NJ: Wiley.

Gao, J., Hu, J., Mao, X., and Perc, M. (2012a). Culturomics meets random fractal theory: insights into long-range correlations of social and natural phenomena over the past two centuries. J. R. Soc. Interface, 7:9. doi: 10.1098/rsif.2011.0846.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gao, J., Hu, J., Mao, X., and Tung, W. (2012b). Detecting low-dimensional chaos by the “noise titration” technique: possible problems and remedies. Chaos Sol. Fract. 45, 213–223. doi: 10.1016/j.chaos.2011.12.004

CrossRef Full Text

Gao, J., Hu, J., and Tung, W. (2011a). Complexity measures of brain wave dynamics. Cogn. Neurodyn. 5, 171–182. doi: 10.1007/s11571-011-9151-3

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gao, J., Hu, J., and Tung, W. (2011b). Facilitating joint chaos and fractal analysis of biosignals through nonlinear adaptive filtering. PLoS ONE 6:e24331. doi: 10.1371/journal.pone.0024331

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gao, J., Hu, J., and Tung, W. (2012c). Entropy measures for biological signal analysis. Nonlin. Dyn. 68, 431–444. doi: 10.1007/s11071-011-0281-2

CrossRef Full Text

Gao, J., Hu, J., Tung, W., and Blasch, E. (2012d). Multiscale analysis of physiological data by scale-dependent lyapunov exponent. Front. Fract. Physiol. 2:110. doi: 10.3389/fphys.2011.00110.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gao, J., Hu, J., Tung, W., and Cao, Y. (2006a). Distinguishing chaos from noise by scale-dependent lyapunov exponent. Phys. Rev. E 74:066204. doi: 10.1103/PhysRevE.74.066204

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gao, J., Hu, J., Tung, W., Cao, Y., Sarshar, N., and Roychowdhury, V. (2006b). Assessment of long range correlation in time series: how to avoid pitfalls. Phys. Rev. E 73:016117. doi: 10.1103/PhysRevE.73.016117

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gao, J., Hu, J., Tung, W., and Zheng, Y. (2011c). Multiscale analysis of economic time series by scale-dependent lyapunov exponent. Quant. Finance 13:2. doi: 10.1080/14697688.2011.580774.

CrossRef Full Text

Gao, J., Sultan, H., Hu, J., and Tung, W. (2010). Denoising nonlinear time series by adaptive filtering and wavelet shrinkage: a comparison. IEEE Signal Proc. Lett. 17, 237–240. doi: 10.1109/LSP.2009.2037773

CrossRef Full Text

Gao, J., Tung, W., and Hu, J. (2009). Quantifying dynamical predictability: the pseudo-ensemble approach (in honor of professor andrew majda's 60th birthday). Chi. Ann. Math. B 30, 569–588. doi: 10.1093/bioinformatics/btl190

CrossRef Full Text

Goldberger, A., and West, B. (1987). Applications of nonlinear dynamics to clinical cardiology. Ann. N.Y. Acad. Sci. 504, 155–212. doi: 10.1111/j.1749-6632.1987.tb48733.x

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ho, K., Moody, G., Peng, C., Mietus, J., Larson, M., Levy, D., et al. (1997). Predicting survival in heart failure cases and controls using fully automated methods for deriving nonlinear and conventional indices of heart rate dynamics. Circulation 96, 842–848. doi: 10.1161/01.CIR.96.3.842

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hollingsworth, K., Hodgson, T., Macgowan, G., Blamire, A., and Newton, J. (2011). Impaired cardiac function in chronic fatigue syndrome measured using magnetic resonance cardiac tagging. J. Inter. Med. 271, 264–270. doi: 10.1111/j.1365-2796.2011.02429.x

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hon, E., and Lee, S. (1963). Electronic evaluations of the fetal heart rate patterns preceding fetal death: further observations. Am. J. Obstet. Gynecol. 87, 814–826.

Pubmed Abstract | Pubmed Full Text

Hu, J., Gao, J., and Tung, W. (2009a). Characterizing heart rate variability by scale-dependent lyapunov exponent. (special issue on Controversial Topics in Nonlinear Science: Is the Normal Heart Rate Chaotic?) Chaos 19:028506. doi: 10.1063/1.3152007

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hu, J., Gao, J., Tung, W., and Cao, Y. (2010). Multiscale analysis of heart rate variability: a comparison of different complexity measures. Ann. Biom. Eng. 38, 854–864. doi: 10.1007/s10439-009-9863-2

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hu, J., Gao, J., and Wang, X. (2009b). Multifractal analysis of sunspot time series: the effects of the 11-year cycle and fourier truncation. J. Stat. Mech. 2009:P02066. doi: 10.1088/1742-5468/2009/02/P02066.

CrossRef Full Text

Hurwitz, B., Coryell, V., Parker, M., Martin, P., Laperriere, A., Klimas, N., et al. (2009). Chronic fatigue syndrome: illness severity, sedentary lifestyle, blood volume and evidence of diminished cardiac function. Clin. Sci. (Lond.) 118, 125–135. doi: 10.1042/CS20090055

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kaplan, D., and Goldberger, A. (1991). Chaos in cardiology. J. Cardiovasc. Electrophysiol. 2, 342–354. doi: 10.1111/j.1540-8167.1991.tb01331.x

CrossRef Full Text

King, D., Ogilvie, M., Pereira, B., Chang, Y., Manning, R., Conner, J., et al. (2009). Heart rate variability as a triage tool in patients with trauma during prehospital helicopter transport. J. Trauma 67, 236–240. doi: 10.1097/TA.0b013e3181ad67de

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kleiger, R., Miller, J., Bigger, J. Jr., and Moss, A. (1987). Decreased heart rate variability and its association with increased mortality after acute myocardial infarction. Am. J. Cardiol. 59, 256–262. doi: 10.1016/0002-9149(87)90795-8

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kobayashi, M., and Musha, T. (1982). 1/f fluctuation of heart beat period. IEEE Trans. Biomed. Eng. 29, 456–457. doi: 10.1529/biophysj.108.139691

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kuznetsov, N., Bonnette, S., Gao, J., and Riley, M. (2012). Adaptive fractal analysis reveals limits to fractal scaling in center of pressure trajectories. Ann. Biomed. Eng. doi: 10.1007/s10439-012-0646-9. [Epub ahead of print].

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Malik, M. (1996). Task force of the european society of cardiology & the north american society of pacing & electrophysiology: heart rate variability: standards of measurement, physiological interpretation, and clinical use. Circulation 93, 1043–1065. doi: 10.1161/01.CIR.93.5.1043

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Miwa, K., and Fujita, M. (2011). Small heart with low cardiac output for orthostatic intolerance in patients with chronic fatigue syndrome. Clin. Cardiol. 34, 782–786. doi: 10.1002/clc.20962

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Naschitz, J., Fields, M., Isseroff, H., Sharif, D., Sabo, E., and Rosner, I. (2006). Shortened qt interval: a distinctive feature of the dysautonomia of chronic fatigue syndrome. J. Electrocardiol. 39, 389–394. doi: 10.1016/j.jelectrocard.2005.10.014

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Newton, J., Pairman, J., Hallsworth, K., Moore, S., Plotz, T., and Trenell, M. (2011). Physical activity intensity but not sedentary activity is reduced in chronic fatigue syndrome and is associated with autonomic regulation. Q. J. Med. 104, 681–687. doi: 10.1093/qjmed/hcr029

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Otzenberger, H., Gronfier, C., Simon, C., Charloux, A., Ehrhart, J., Piquard, F., et al. (1998). Dynamic heart rate variability: a tool for exploring sympathovagal balance continuously during sleep in men. Am. J. Physiol. 275, H946–H950.

Pubmed Abstract | Pubmed Full Text

Packard, N., Crutchfield, J., Farmer, J., and Shaw, R. (1980). Gemomtry from time-series. Phys. Rev. Lett. 45, 712–716. doi: 10.1103/PhysRevLett.45.712

CrossRef Full Text

Peng, C., Buldyrev, S., Havlin, S., Simons, M., Stanley, H., and Goldberger, A. (1994). On the mo-saic organization of dna sequences. Phys. Rev. E 49, 1685–1689. doi: 10.1103/PhysRevE.49.1685

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Peng, C., Mietus, J., Hausdorff, J., Havlin, S., Stanley, H., and Goldberger, A. (1993). Long-range anticorrelations and non-gaussian behavior of the heartbeat. Phys. Rev. Lett. 70, 1343–1346. doi: 10.1103/PhysRevLett.70.1343

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Pincus, S., and Viscarello, R. (1992). Approximate entropy: a regularity statistic for fetal heart rate analysis. Obst. Gynecol. 79, 249–255.

Pubmed Abstract | Pubmed Full Text

Rahman, K., Burton, A., Galbraith, S., Lloyd, A., and Vollmer-Conna, U. (2011). Sleep-wake behavior in chronic fatigue syndrome. SLEEP 34, 671–678.

Pubmed Abstract | Pubmed Full Text

Reeves, W., Wagner, D., Nisenbaum, R., Jones, J., Gurbaxani, B., Solomon, L., et al. (2005). Chronic fatigue syndrome: a clinically empirical approach to its definition and study. BMC Med. 3:19 doi: 10.1186/1741-7015-3-19

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Riley, M., Kuznetsov, N., Bonnette, S., Wallot, W., and Gao, J. (2012). A tutorial introduction to adaptive fractal analysis. Front. Fract. Physiol. 3:371. doi: 10.3389/fphys.2012.00371.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ryan, D., and Sarson, G. (2008). The geodynamo as a low-dimensional deterministic system at the edge of chaos. Europhys. Lett. 83:49001. doi: 10.1209/0295-5075/83/49001

CrossRef Full Text

Sauer, T., Yorke, J., and Casdagli, M. (1991). Embedology. J. Stat. Phys. 65, 579–616. doi: 10.1007/BF01053745

CrossRef Full Text

Smith, A.-L., and Reynolds, K. (2006). Survey of poincare indices for measuring heart rate variability. Aust. Phys. Eng. Sci. Med. 29, 97–101.

Pubmed Abstract | Pubmed Full Text

Smith, A.-L., Reynolds, K., and Owen, H. (2007). Correlated poincare indices for measuring heart rate variability. Aust. Phys. Eng. Sci. Med. 30, 336–341.

Pubmed Abstract | Pubmed Full Text

Takens, F. (1981). “Detecting strange attractors in turbulence,” in Dynamical Systems and Turbulence, Lecture Notes in Mathematics. Vol. 898, eds D. A. Rand, L. S. Young (Berlin: Springer-Verlag), 366.

Tung, W., Gao, J., Hu, J., and Yang, L. (2011). Recovering chaotic signals in heavy noise environments. Phys. Rev. E 83:046210. doi: 10.1103/PhysRevE.83.046210

CrossRef Full Text

West, B. (2010). Fractal physiology and the fractional calculus: a perspective. Front. Physiol. 1, 1–17. doi: 10.3389/fphys.2010.00012

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text