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).

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.

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⁽ⁱ⁺¹⁾(l₂), respectively, where l₁, l₂ = 1, …, 2n + 1. We then define the fitting for the overlapped region as (1)y(c)(l)=w1y(i)(l+n)+w2y(i+1)(l), l=1, 2, …, n+1, where w1=(1−l−1n) and w2=l−1n 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) (2)F(w)=[1N∑i=1N(u(i)−v(i))2]1/2~wH. 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.

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: (3)Vi=[x(i), x(i+L), …, x(i+(m−1)L)],  i=1, …,      Np=n−(m−1)L, 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, (4)εt+Δt=εteλ(εt)Δt, where λ (ε_t) is the SDLE given by (5)λ(εt)=lnεt+Δt−lnεtΔt.

Equivalently, we can express this as, (6)dεtdt=λ(εt)εt.

To compute SDLE, we check whether pairs of vectors (V_i, V_i) defined by Equation (3) satisfy the following Inequality, (7)εk≤‖Vi−Vj‖≤εk+Δεk,   k=1, 2, 3, …, where ε_k and Δε_k are arbitrarily chosen small distances, and (8)||Vi−Vj||=∑w=1m(xi+(w−1)L−xj+(w−1)L)2

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 (9)λ(εt)= 〈ln‖Vi+t+Δt−Vj+t+Δt‖−ln‖Vi+t−Vj+t‖〉Δt 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 (10)εt=‖Vi+t−Vj+t‖=∑w=1m(xi+(w−1)L+t−xj+(w−1)L+t)2

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

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.

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.

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.

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:

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).

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).

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.