Using local scale exponent to characterize heart rate variability in response to postural changes in people with spinal cord injury

Heart rate variability (HRV) is a promising marker for evaluating the remaining autonomic function in people with spinal cord injury (SCI). HRV is commonly assessed by spectral analysis and detrended fluctuation analysis (DFA). This study aimed to investigate whether local scale exponent α(t) can reveal new features of HRV that cannot be reflected by spectral measures and DFA coefficients. We studied 12 participants with SCI and 15 healthy able-bodied controls. ECG signals were continually recorded during 10 min sitting and 10 min prone postures. α(t) was calculated for scales between 4 and 60 s. Because α(t) could be overestimated at small scales, we developed an approach for correcting α(t) based on previous studies. The simulation results on simulated monofractal time series with α between 0.5 and 1.3 showed that the proposed method can yield improved estimation of α(t). We applied the proposed method to raw RR interval series. The results showed that α(t) in healthy controls monotonically decreased with scale at scales between 4 and 12 s (0.083–0.25 Hz) in both the sitting and prone postures, whereas in participants with SCI, α(t) slowly decreased at almost all scales. The sharp decreasing trend in α(t) in controls suggests a more complex dynamics of HRV in controls. α(t) at scales between 4 (0.25 Hz) and around 7 s (0.143 Hz) was lower in subjects with SCI than in controls in the sitting posture; α(t) at a narrow range of scales around 12 s (0.083 Hz) was higher in participants with SCI than in controls in the prone posture. However, none of normalized low frequency (0.04–0.15 Hz) power, the ratio of low frequency power to high frequency (0.15–0.4 Hz) power and long-term (>11 beats) DFA coefficient showed significant difference between healthy controls and subjects with SCI in the prone posture. Our results suggest that α(t) can reveal more detailed information in comparison to spectral measures and the standard DFA parameters.

Introduction

Spinal cord injury (SCI) disrupts the descending autonomic pathways, resulting in a variety of dysfunction of the cardiovascular system related to the level and severity of injury to these pathways (Alexander et al., 2009). This commonly manifests as neurogenic shock, orthostatic hypotension, autonomic dysreflexia, and cardiac rhythm disturbances (Claydon and Krassioukov, 2006; Krassioukov and Claydon, 2006). To date, the assessment to quantify severity of injury to autonomic pathways has not been well developed, and there is still no consensus on this issue (West et al., 2013). The clinical practice focuses on the assessment of motor and sensory impairments following SCI according to the International Standards for the Neurological Classification of SCI, rather than assessment of severity of injury to autonomic pathways (Claydon and Krassioukov, 2008; Jan et al., 2013).

Heart rate variability (HRV) represents one of the most promising markers of autonomic activities (Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology, 1996). It has been widely accepted that fluctuations in heart rate (HR) occurring at different frequencies reflect the activities of autonomic neural outflows (West et al., 2013). Previous research has shown great potential of HRV for evaluating the remaining autonomic functions of the cardiovascular system in patients with SCI (Merati et al., 2006; Claydon and Krassioukov, 2008; Jan et al., 2013). Various methods have been developed to quantify HRV, including time domain, frequency domain, and nonlinear methods (Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology, 1996). The power spectra density (PSD) of HRV reveals two characteristic frequencies: one is defined as the low frequency (LF, 0.04–0.15 Hz) and the other is defined as the high frequency (HF, 0.15–0.4 Hz). HF of HRV reflects vagal modulation, and LF is considered to be jointly mediated by sympathetic and vagal nerves (Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology, 1996). The ratio of LF to HF power is widely used as an index of sympathovagal balance for assessing cardiovascular regulation. Indeed, time and frequency domain indexes of HRV have provided prognostic information in many pathological conditions. However, in addition to autonomic modulations, HRV can also be influenced by other mechanisms, e.g., intrinsic variations of pacemaker rate (Ponard et al., 2007) or fluctuations in circulating neurohumoral factors (Galetta et al., 2008). These mechanisms may operate at time scales that time and frequency domain methods cannot capture adequately (Tan et al., 2009). Hence, scale independent measures have been introduced to study the integrated control of HR. The use of such measure is based on the observation that HR fluctuations exhibit scale-invariant patterns over a wide range of time scales that breakdown under pathological conditions (Peng et al., 1995; Goldberger et al., 2002).

Detrended fluctuation analysis (DFA) (Peng et al., 1995) is one of the most commonly used scale independent method. DFA provides a quantitative parameter, the scaling exponent α, to represent the correlation properties of the data (Peng et al., 1995). Usually, the properties of HRV data (RR interval series) are described by two scaling exponents: a short-term (4–11 beats) exponent α₁ and a long-term (>11 beats) exponent α₂ (Beckers et al., 2006). It was demonstrated that α₁ of HRV can reflect changes in autonomic tone induced by exercise, maneuvers such as passive head-up tilt, cold hand immersion, cold face immersion (Tulppo et al., 2001, 2005), and aging or cardiac pathologies (Beckers et al., 2006). A growing number of studies have suggested that α₁ can yield prognostic information that cannot be provided by conventional measures (Huikuri et al., 2009). However, the use of a single scaling exponent for characterizing the HR dynamics has been questioned (Ivanov et al., 1999; Echeverria et al., 2003; Castiglioni et al., 2009, 2011). It has been recognized that HRV does not always present a uniform power-law behavior, especially in abnormal physiological conditions (Echeverria et al., 2003). Ivanov et al. (1999) suggested that human heartbeat dynamics requires a large number of exponents to characterize its properties. A recent study even suggested that the physiological effects of autonomic outflow may mask intrinsic fractal behavior of the sinoatrial node (Tan et al., 2009). Some authors therefore proposed to use a whole spectrum of local scale exponents α(t) to describe HR dynamics (Bojorges-Valdez et al., 2007; Castiglioni et al., 2009, 2011), where t is the time scale. Castiglioni et al. (2011) suggested that α(t) describes the local correlation properties of data. By using autonomic blocking agents, they showed that both cardiac sympathetic and vagal outflows affect the spectrum of α(t) but with opposite effects. However, it is still unclear whether the spectrum of α(t) reveals new features that are not characterized by the PSD and DFA coefficients. In addition, α(t) at small scales can be overestimated (Viswanathan et al., 1997; Kantelhardt et al., 2001). This problem was not well addressed in the previous studies (Castiglioni et al., 2009, 2011). This is an important issue because α(t) at small scales and their mean value (i.e., α₁), have been used to assess HR dynamics due to their sensitivity to changes in cardiac autonomic tone (Tulppo et al., 2001, 2005).

The purpose of this study was to investigate whether α(t) reveals new features of HRV in people with SCI that cannot be reflected by spectral measure and DFA coefficients. This was done by changing the posture from the sitting to prone posture in people with SCI, which is known to induce autonomic modulated vasodilatory responses. To accurately estimate α(t), we developed an approach for correcting α(t) at small scales based on a previous study by Kantelhardt et al. (2001). Our method may be used to evaluate the local properties of HRV and may enhance the current understanding of residual autonomic function in people with SCI.

Methods

Subjects

Twenty seven participants were recruited in this study, including 12 people with SCI and 15 healthy able-bodied controls. Their demographic data are shown in Table 1. The SCI group included five participants with a spinal injury level at C4–T5 [five incomplete, American Spinal Injury Association Impairment Scale (AIS) B, C, or D] and seven subjects with a spinal injury level at T6–T12 [four complete (AIS A) and three incomplete (AIS B, C, or D)]. All participants with SCI were in a stable clinical condition (the injury event occurred more than 6 months before the time of the study). None of the subjects had any diagnosed cardiovascular or neurological diseases that might affect autonomic cardiovascular control. All participants gave informed consent to participate in this study, which was approved by the Institutional Review Board of the University of Oklahoma Health Sciences Center.

TABLE 1

Table 1. Demographic data of the enrolled subjects.

Data Collection

All experiments were performed in a university laboratory. Room temperature was maintained at about 23°C. Each subject was asked to stay in the laboratory for at least 30 min to acclimate to the room temperature. When a subject sat in the wheelchair, the electrocardiography (ECG) electrodes of a three-lead Biopac ECG monitor (ECG100C, Biopac Systems; Goleta, California) were placed on the right ventral wrist, right medial ankle, and left medial ankle. The ECG signals were recorded for 10 min with a sampling rate of 1000 Hz. Then the subject was transferred to a mat table in a prone position for 10 min recording. Able-bodied controls followed the same procedures but changed from the sitting posture to the prone posture without assistance. By using the AcqKnowledge software (Biopac Systems), the ECG signal was edited to remove artifacts and then filtered with a bandpass filter of 0.5–32 Hz. After RR peak detection and visual inspection by the operator, the consecutive RR intervals (RRIs) were exported for later processing.

Linear Analysis

Mean and standard deviation (SD) of RRI series were calculated in the time domain. Then the RRI series was resampled at 2 Hz using a cubic spline approximation (Aubert et al., 1999). The PSD of HRV was calculated using an autoregressive model with afixed order of 16 (Boardman et al., 2002). Three parameters of PSD were calculated according to the published guideline (Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology, 1996): normalized low frequency power [% LF = 100×LF power/(total power-VLF power)], high frequency power, and LF to HF ratio (LF power/HF power), where total power was defined as the variance of the data series and VLF was the very low frequency range (<0.04 Hz).

Detrended Fluctuation Analysis (DFA)

The DFA method applied to RRI series has been described elsewhere (Peng et al., 1995). We have briefly described the methods used in this study here. Given a series of N RRIs, {x(i), i = 1,…,N}, it was first integrated after mean subtraction

$\begin{array}{cc}y\left(k\right)={\displaystyle {\sum }_{i=1}^k\left(x\left(i\right)-\langle RR\rangle \right)},& (1)\end{array}$

where 〈RR〉 was the mean of the series. Next, the integrated series y(k) was divided into boxes of n RRIs. In each box, the local trend, y_n(k), was estimated by a least-square fit of y(k) using a polynomial function (a linear function was used throughout this work). The root-mean square fluctuation

$\begin{array}{cc}F\left(n\right)=\sqrt{\frac{1}{N}{\displaystyle {\sum }_{k=1}^N({(y(k)-y_n\left(k\right))}^2)}}& (2)\end{array}$

was calculated for box sizes n≥ 4. A power-law between F(n) and n indicated the presence of scaling: F(n) ~ n^α. The parameter α, called scaling exponent, was estimated by the slope of the log[F(n)]–log(n) plot. α represents the correlation properties of the signal. For examples, α = 0.5 indicates white noise, α = 1 represents the behavior of a 1/f process having persistent long range correlations, and α = 1.5 indicates Brownian noise. For HRV data, standard DFA typically includes a short-term (4–11 beats) scaling exponent α₁ and a long-term (>11 beats) scaling exponent α₂ (Beckers et al., 2006).

Local Scale Exponent

A local scale exponent was defined as the derivative of log[F(n)] with respect to log(n)

$\begin{array}{cc}\alpha (n_k)=\frac{log\left[F(n_{k+1})\right]-log[F(n_{k-1})]}{log(n_{k+1})-log(n_{k-1})},& (3)\end{array}$

where {n_k} was a set of box sizes spaced evenly on a log scale. An intrinsic problem was that the log[F(n)] − log(n) plot deviates from scaling at small scales, resulting in over-estimations of scale exponents at small scales (Viswanathan et al., 1997; Kantelhardt et al., 2001; Echeverria et al., 2003). To address this problem, Kantelhardt et al. (2001) introduced a correction function

$\begin{array}{cc}{K}^{\alpha }(n)=\frac{{\langle {[F(n)]}^2\rangle }^{1/2}/n^{\alpha }}{{\langle {[F(n\prime )]}^2\rangle }^{1/2}/n{\prime }^{\alpha }},& (4)\end{array}$

where 〈·〉 denotes the average over different configurations, and n′ is a specific box size. Because K^α(n) depends only weakly on α, K^α(n) for uncorrelated data, i.e., K^1/2(n), can be used in all cases. K^1/2(n) can be obtained by analyzing the corresponding shuffled data, where all long-range correlations are destroyed by shuffling. Kantelhardt et al. (2001) suggested that n′ has to be large (n′ > 50) but must remain significantly smaller than the record length N and they suggested n′ ≈ N/20 to be a reasonable number. The modified fluctuation function is given by Kantelhardt et al. (2001)

$\begin{array}{cc}{F}_{mod}\left(n\right)=F(n)/K^{1/2}(n).& (5)\end{array}$

However, we found that for short time series, this method could result in underestimated α(n_k) at small scales even for using a linear function y_(k) when calculating the fluctuation function F(n) (Figure 1). Furthermore, it seems that a larger value of α leads to more prominent bias. The reason is that a larger value of α leads to more prominent differences between K^α(n) (Equation 4) and K^1/2(n) at small scales (see Figure 2). Thus, the underestimation bias in α(n) can be corrected as

$\begin{array}{cc}\alpha \left(n\right)={\alpha }_K(n)·K^{\alpha }(n)/K^{\frac{1}{2}}(n),& (6)\end{array}$

where α_K(n) denotes local scale exponents obtained by using the method by Kantelhardt et al.

FIGURE 1

Figure 1. Local scale exponents α(n) for monofractal time series with α between 0.5 and 1.3. The monofractal time series were generated by inverse Fourier transform of the power spectrum S(f) ∝ f^−β with β = 2α −1. When calculating α(n), linear trends were removed. The results were obtained by averaging over 50 series having a length of 750 data points. The circles represent the values without correction and the stars represent the values corrected using the method proposed by Kantelhardt et al. (2001). The error bars represent the standard deviations.

FIGURE 2

Figure 2. The correction function K^α(n) for α = 0.5, 0.7, and 1.3. A larger value of α leads to more prominent deviation between K^α(n) (Equation 4) and K^1/2(n). The error bars represent the standard deviations.

To verify whether the proposed method improves the estimation of α(n), we considered simulated monofractal time series with α between 0.5 and 1.3, which were generated by inverse Fourier transform of the power spectrum S(f) ∝ f^−β with β = 2α −1. In Figure 3, the values of α(n) obtained by using Equation (6) are compared to those obtained by using Equation (5). Obviously, the suggested method substantially improves estimations of α(n) at small scales. Since for RRI series, the mean value of α(n) at n < 11 beats, i.e., the short-term DFA coefficient α₁, is generally larger than 1 (Beckers et al., 2006), we argue that the proposed method can yield improved results.

FIGURE 3

Figure 3. Local scale exponent α(n) for monofractal time series obtained by using the method proposed by Kantelhardt et al. (Kantelhardt correction) and by the proposed method. In these calculations, linear trends were removed. The results were obtained by averaging over 50 series having a length of 750 data points. The error bars represent the standard deviations.

For a RR interval series, the procedure of the proposed correction method is as follows. (1) Calculation of the short-term DFA coefficient α1 of the RR interval series. In general, the scales for calculating α₁ range between four beats and a specific scale around 11 beats. (2) Estimation of the actual scaling exponent α from α₁. Figure 4 shows the relationship between the actual exponent α and α₁ for monofractal time series with α between 0.5 and 1.5, where the monofractal time series of the same length of the RR interval series can be generated by inverse Fourier transform. According to this relationship, i.e., the fitting curve (in red color), it is easy to estimate α from α₁. (3) Calculation of local scale exponents using the method proposed by Kantelhardt et al. (4) Calculation of local scale exponent using Equation (6).

FIGURE 4

Figure 4. Relationship between actual scaling exponent α and DFA coefficient at scales between 4 and 11. The circles represent the mean values of DFA coefficients of 50 simulated monofractal signals with a length of 750 points. The lengths of the error bars represent the standard deviations. The fitting curve (red line) was computed in a least-squares sense.

Application to HRV Data

Because the time duration of a box size depends on the mean RRI, if two series have different values of mean RRI, the same box size actually represents different time scales. This is the case in many studies under different experimental conditions. In this study, mean RRI was generally shorter in the sitting posture compared to the prone posture, especially in able-bodied controls. In order to compare scale exponent spectrum of HRV between the sitting and prone postures and among different subjects, we converted the box sizes {n_k} to time scales {t_k} by Castiglioni et al. (2009)

$\begin{array}{cc}{t}_k=n_k\times \langle RR\rangle .& (7)\end{array}$

Then the spectrum α(t_k) was linearly interpolated to obtain a new spectrum α(t_i) with t_i evenly distributed between 4 and 60 s on a log scale.

The rationales for choosing the time scale interval are as follows. In DFA of RRI series, the smallest box size is usually chosen as four beats (Beckers et al., 2006). For most of our data sets, mean RRI ranges between 0.6 and 1 s. Thus, the time scale corresponding to the smallest box size of four RRIs ranges between 2.4 and 4 s. We therefore chose 4 s as the smallest time scale. In the previous studies by Castiglioni et al. (2009, 2011), the smallest time scale was also chosen as 4 s. On the other hand, the largest time scale depends on the length of data series (10 min). Kantelhardt (2009) suggested that the largest scale can be N/4. Castiglioni et al. (2011) reported that α(n) can be estimated without substantial bias up to scales equal to N/7. In their study, the largest time scale for 15–20 min record was chosen as 100 s. However, our simulation results on monofractal time series indicated that, when scales exceed N/20, α(n) may become unstable and deviate from the expected values (Figures 1, 3). Thus, we speculate that for RRI series derived from 10 min ECG signals, the reliability of estimations of α(n) decreases when time scale exceeds 30 s. The time scale range 4–30 s corresponds to the frequency interval 0.033–0.25 Hz. It has been well known that in spectral analysis of HRV, the low frequency (LF) and high frequency (HF) bands are defined as 0.04–0.15 Hz and 0.15–0.4 Hz, respectively. According to the published guidelines (Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology, 1996), for short-term recordings (~5 min), assessment of very low frequency, i.e., <0.04 Hz, should be avoided. In this study, we also presented α at scales between 30 and 60 s just for reference.

Statistical Analysis

The Wilcoxon signed-ranked test was used to examine the differences of linear measures of HRV, scaling exponents α₁and α₂, and local scaling exponent α(t) between sitting and prone postures (within-subjects test). The Wilcoxon rank-sum test was used to examine the differences in these measures between patients with SCI and able-bodied controls (between-subjects test).

Results

Linear Parameters

Figure 5 shows typical examples of changes in RRI and PSD of HRV in response to the postural change in a healthy able-bodied control (Figures 5A,B,E) and a person with SCI (Figures 5C,D,F). In the control, RRI showed a distinctive increase after the change from the sitting to prone posture (Figures 5A,B). However, RRI in the person with SCI showed only a slight increase after the postural change (Figures 5C,B). The PSD of HRV in the control revealed a decrease in LF power and an increase in HF power (Figure 5E), whereas the PSD of HRV in the person with SCI showed only slight changes (Figure 5F).

FIGURE 5

Figure 5. Typical examples of changes in RRI and PSD of HRV in response to the postural change in an able-bodied control (A,B,E) and a patient with SCI (C,D,F).

Table 2 compares the linear parameters between healthy able-bodied controls and subjects with SCI. Mean RRI was lower in the sitting position in both groups. In controls, %LF was significantly higher, HF power was significantly lower, and thus LF/HF was significantly higher in the sitting posture, whereas in subjects with SCI, these parameters did not show significant change after the postural change. Compared to controls, %LF was significantly lower, HF power was significantly higher, and LF/HF was significantly lower in subjects with SCI in the sitting position. No significant difference was observed between two groups in the prone position.

TABLE 2

Table 2. Linear parameters of HRV in healthy able-bodied controls and subjects with SCI.

DFA Coefficients

The short-term scaling exponent α₁was significantly higher in the sitting posture than in the prone posture in healthy able-bodied controls (p < 0.05) but not in subjects with SCI (Figure 6). In the sitting posture, α₁ was significantly higher in controls than in subjects with SCI (Figure 6A). Unlike α₁, the long-term scaling exponent α₂ showed only slight changes after the postural change in either control or SCI group (Figure 6B) and no significant differences were observed between two groups (in the sitting posture, p = 0.12; in the prone posture, p = 0.10).

FIGURE 6

Figure 6. DFA coefficients α₁ (short-term) and α₂ (long-term) in able-bodied controls and subjects with SCI. Values are means ± SE. *p < 0.05 for within-subjects test; ⁺p < 0.05 for between-subjects test. (A) α₁ in control and SCI groups. (B) α₂ in control and SCI groups.

Local Scale Exponent

Figure 7 presents the local scale exponents α(t) over healthy able-bodied control and SCI groups as means ±SE. In controls in both the sitting and prone postures, α(t) monotonically decreased with scale at the scales between 4 and 12 s and then showed a plateau (Figure 7A). At the scales shorter than 10 s, α(t) was significantly higher in the sitting posture than in the prone posture (Figure 7A), while at the scales larger than 15 s, α(t) remained at lower values for both postures. In subjects with SCI, α(t) slowly decreased at almost all scales. No differences were observed between two postures (Figure 7B). When comparing controls and subjects with SCI, in the sitting posture, α(t) in subjects with SCI was significantly lower at scales shorter than 8 s and tended to be higher at larger scales (Figure 7C). In the prone posture, α(t) at a narrow range of scales around 12 s was significantly higher in subjects with SCI than in controls (Figure 7D).

FIGURE 7

Figure 7. Local scale exponents αt in healthy able-bodied controls and subjects with SCI. Values are presented as means ± SE. *p < 0.05 for within-subjects test; ⁺p < 0.05 for between-subjects test. (A) Comparisons of α(t) in control group between two postures. (B) Comparisons of α(t) in SCI group between two postures. (C) Comparison of α(t) in the sitting position between two groups. (D) Comparisons of α(t) in the prone position between two groups.

Discussion

The main findings of this study are: (1) local scale exponent α(t) in healthy able-bodied controls rapidly decreased with scale at small scales in both the sitting and prone postures but slowly decreased in subjects with SCI (Figures 7A,B); (2) in the sitting posture, α(t) at small scales was lower in subjects with SCI than in controls (Figure 7C); and (3) in the prone posture, α(t) at moderate scales was higher in subjects with SCI than in controls. However, normalized low frequency (0.04–0.15 Hz) power, the ratio of low frequency power to high frequency (0.15–0.4 Hz) power and long-term (>11 beats) DFA coefficient α₂ did not show significant difference between healthy controls and subjects with SCI in the prone posture (Table 2, Figure 6). Our findings support our hypothesis that α(t) can reveal important features of HRV in patients with SCI that are not reflected by the sympathovagal balance and DFA coefficients. This approach may be used to investigate the effects of SCI-induced autonomic damage on cardiovascular dysfunction.

Spectral analysis of HRV is considered to be a useful tool for evaluating cardiovascular autonomic control after SCI (Claydon and Krassioukov, 2008; Jan et al., 2013). In the present study, sympathovagal balance in controls increased in the sitting posture and decreased in the prone posture, whereas in patients with SCI, it showed only small changes (Table 2). These results were consistent with the previous studies (Claydon and Krassioukov, 2008; Jan et al., 2013). Increased sympathovagal balance in the sitting posture is associated with withdraw of vagal activity and enhanced sympathetic outflow to the heart (Wecht et al., 2006; Claydon and Krassioukov, 2008). In subjects with SCI, the increase in sympathovagal balance was attenuated, probably reflecting limited ability of vagal withdrawal to influence HR and reduced sympathetic cardiac modulation (Wecht et al., 2006). On the other hand, in the sitting posture, the short-term DFA coefficient α₁ significantly increased in controls but not in subjects with SCI (Figure 6). These results were similar to those reported by Tulppo et al. (2001) and Merati et al. (2006). Our results suggested that the LF to HF because it has been demonstrated that ratio and DFA coefficients provide similar characterizations of HRV. This is not surprising, DFA coefficients can be obtained from the PSD (Willson et al., 2002). Willson et al. (2002) suggested that α₁ is related to 2/[1 + (HF/LF)] and α₂ is related to 2/[1 + (LF/VLF)] with VLF representing the power in the very low frequency region below 0.04 Hz. Platisa and Gal (2006) observed an approximated linear relationship between α₁ and ln(LF/HF). Thus, it seems that α₁ reflects the balance between the LF and HF oscillations and α₂ reflects the balance between the LF and VLF oscillations.

Unlike DFA coefficients, local scale exponent reflects local properties of HRV. Castiglioni et al. (2011) suggested that local scale exponent reflects the local correlation properties of the data. By using autonomic blocking agents, they showed that the vagal outflow contributes with white-noise components to HRV while the cardiac sympathetic outflow adds Brownian motion components at short scales and contributes to a plateau between 40 and 80 s. In the present study, in able-bodied controls, α(t) at t between 4 and 12 s was significantly higher in the sitting posture than in the prone posture (Figure 7A). These scales correspond to frequencies between 0.083 and 0.25 Hz, which spread across the LF and HF bands. Supposing that the mean RRI is 0.8 s (see Table 2), the scales between 4 and 12 s correspond to 5–15 beats (Eq. 7). This range of beat number may spread across the ranges of observation window sizes over which α₁ and α₂ are calculated. As shown in Figure 6, α₂ in controls in the sitting posture tended to be higher than in the prone posture but the differences did not reach a significant level. The reason might be that α₂ was calculated over a wide range of observation window sizes, while in major parts of the range, α(t) exhibited similar behavior under two conditions (Figure 7A).

In the sitting posture, α(t) at t < 8 s was significantly lower in subjects with SCI than in controls (Figure 7C). The lower values of α(t) are compatible with lower values of α₁ (Figure 6A) and lower values of sympathovagal balance (Table 2). At t between 6 and 15 s, roughly corresponding to 0.07–0.16 Hz, α(t) showed different curve configurations between two groups: it sharply decreased with t in controls but slowly decreased in subjects with SCI. The sharp decreasing trend in α(t) in controls suggests a more complex dynamics of HRV in controls. In contrast, the relative constant values of α(t) in subjects with SCI might reflect an attenuation of HRV. These features were not reflected by sympathovagal balance quantified by the spectral analysis and DFA coefficients.

In the prone posture, α(t) at t around 12 s was significantly higher in subjects with SCI than in controls (Figure 7D). The higher values of α(t) are compatible with the higher values of α₂ in subjects with SCI although the differences in α₂ between two groups did not reach the significant level (Figure 6B). The reasons for the higher values of α(t) and α₂ in subjects with SCI are not clear. A possible explanation is the various levels and severities of SCI. It is well known that the mechanisms of HR regulation after SCI depend on the level and severity of injury (Claydon and Krassioukov, 2008). For example, it was reported that in the supine position, sympathovagal balance in controls was higher than in subjects with high-level lesions but lower than in those with low-level lesions (Merati et al., 2006; Claydon and Krassioukov, 2008). Accordingly, it was suggested that at rest, reduced sympathetic outflow are not balanced by reduced vagal outflow in subjects with cervical SCI and mechanisms underlying reduced vagal tone in subjects with thoracic SCI subjects are uncertain (Claydon and Krassioukov, 2008). Nevertheless, a striking phenomenon was that α(t) maintained relatively constant in subjects with SCI, whereas it monotonically decreased with scale in controls (Figure 7D). This might suggest that HRV in subjects with SCI was less complex compared to controls.

This study has several limitations. First, we only recruited 12 subjects with SCI into this study. However, our results confirmed our hypothesis that the spectrum of α(t) reveals important features of HRV in subjects with SCI that are not reflected by the spectral analysis (LF to HF ratio) and DFA coefficients. Second, we did not investigate the influences of injury level on α(t), because the C4–T5 group mainly consisted of incomplete spinal injury, while the T6–T12 group mainly consisted of complete spinal injury. Further studies may be required to determine the influences of injury level and completeness on α(t). Third, the SCI group was older compared to the control group. The age issue may affect the validity of our results. However, the age effect is a minor confounding variable when compared to postural changes (Jan et al., 2013). Further studies should consider using an age-matched research design.

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

This study was supported in part by the Paralyzed Veterans of America Research Foundation (PVA2827).

References

Alexander, M. S., Biering-Sorensen, F., Bodner, D., Brackett, N. L., Cardenas, D., Charlifue, S., et al. (2009). International standards to document remaining autonomic function after spinal cord injury. Spinal Cord 47, 36–43. doi: 10.1038/sc.2008.121

PubMed Abstract | CrossRef Full Text | Google Scholar

Aubert, A. E., Ramaekers, D., Beckers, F., Breem, R., Denef, C., Van de Werf, F., et al. (1999). The analysis of heart rate variability in unrestrained rats. Validation of method and results. Comput. Methods Programs Biomed. 60, 197–213. doi: 10.1016/S0169-2607(99)00017-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Beckers, F., Verheyden, B., and Aubert, A. E. (2006). Aging and nonlinear heart rate control in a healthy population. Am. J. Physiol. Heart Circ. Physiol. 290, H2560–H2570. doi: 10.1152/ajpheart.00903.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Boardman, A., Schlindwein, F. S., Rocha, A. P., and Leite, A. (2002). A study on the optimum order of autoregressive models for heart rate variability. Physiol. Meas. 23, 325–336. doi: 10.1088/0967-3334/23/2/308

PubMed Abstract | CrossRef Full Text | Google Scholar

Bojorges-Valdez, E. R., Echeverria, J. C., Valdes-Cristerna, R., and Pena, M. A. (2007). Scaling patterns of heart rate variability data. Physiol. Meas. 28, 721–730. doi: 10.1088/0967-3334/28/6/010

PubMed Abstract | CrossRef Full Text | Google Scholar

Castiglioni, P., Parati, G., Civijian, A., Quintin, L., and Di Rienzo, M. (2009). Local scale exponents of blood pressure and heart rate variability by detrended fluctuation analysis: effects of posture, exercise, and aging. IEEE Trans. Biomed. Eng. 56, 675–684. doi: 10.1109/TBME.2008.2005949

PubMed Abstract | CrossRef Full Text | Google Scholar

Castiglioni, P., Parati, G., Di Rienzo, M., Carabalona, R., Cividjian, A., and Quintin, L. (2011). Scale exponents of blood pressure and heart rate during autonomic blockade as assessed by detrended fluctuation analysis. J. Physiol. 589(Pt 2), 355–369. doi: 10.1113/jphysiol.2010.196428

PubMed Abstract | CrossRef Full Text | Google Scholar

Claydon, V. E., and Krassioukov, A. V. (2006). Orthostatic hypotension and autonomic pathways after spinal cord injury. J. Neurotrauma 23, 1713–1725. doi: 10.1089/neu.2006.23.1713

PubMed Abstract | CrossRef Full Text | Google Scholar

Claydon, V. E., and Krassioukov, A. V. (2008). Clinical correlates of frequency analyses of cardiovascular control after spinal cord injury. Am. J. Physiol. Heart Circ. Physiol. 294, H668–H678. doi: 10.1152/ajpheart.00869.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Echeverria, J. C., Woolfson, M. S., Crowe, J. A., Hayes-Gill, B. R., Croaker, G. D., and Vyas, H. (2003). Interpretation of heart rate variability via detrended fluctuation analysis and alphabeta filter. Chaos 13, 467–475. doi: 10.1063/1.1562051

PubMed Abstract | CrossRef Full Text | Google Scholar

Galetta, F., Franzoni, F., Fallahi, P., Tocchini, L., Braccini, L., Santoro, G., et al. (2008). Changes in heart rate variability and QT dispersion in patients with overt hypothyroidism. Eur. J. Endocrinol. 158, 85–90. doi: 10.1530/EJE-07-0357

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldberger, A. L., Peng, C. K., and Lipsitz, L. A. (2002). What is physiologic complexity and how does it change with aging and disease? Neurobiol. Aging 23, 23–26. doi: 10.1016/S0197-4580(01)00266-4

PubMed Abstract | CrossRef Full Text

Huikuri, H. V., Perkiomaki, J. S., Maestri, R., and Pinna, G. D. (2009). Clinical impact of evaluation of cardiovascular control by novel methods of heart rate dynamics. Philos. Trans. A Math. Phys. Eng. Sci. 367, 1223–1238. doi: 10.1098/rsta.2008.0294

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanov, P. C., Amaral, L. A., Goldberger, A. L., Havlin, S., Rosenblum, M. G., Struzik, Z. R., et al. (1999). Multifractality in human heartbeat dynamics. Nature 399, 461–465. doi: 10.1038/20924

PubMed Abstract | CrossRef Full Text | Google Scholar

Jan, Y. K., Anderson, M., Soltani, J., Burns, S., and Foreman, R. D. (2013). Comparison of changes in heart rate variability and sacral skin perfusion in response to postural changes in people with spinal cord injury. J. Rehabil. Res. Dev. 50, 203–214. doi: 10.1682/JRRD.2011.08.0138

PubMed Abstract | CrossRef Full Text | Google Scholar

Kantelhardt, J. W. (2009). “Fractal and multifractal time series,” in Encyclopedia of Complexity and System Science, ed R. A. Meyers (New York, NY: Springer), 3754–3779.

Google Scholar

Kantelhardt, J. W., Koscielny-Bunde, E., Rego, H. A., Havlin, S., and Bunde, A. (2001). Detecting long-range correlations with detrended fluctuation analysis. Phys. A 295, 441–454. doi: 10.1016/S0378-4371(01)00144-3

CrossRef Full Text | Google Scholar

Krassioukov, A., and Claydon, V. E. (2006). The clinical problems in cardiovascular control following spinal cord injury: an overview. Prog. Brain Res. 152, 223–229. doi: 10.1016/S0079-6123(05)52014-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Merati, G., Di Rienzo, M., Parati, G., Veicsteinas, A., and Castiglioni, P. (2006). Assessment of the autonomic control of heart rate variability in healthy and spinal-cord injured subjects: contribution of different complexity-based estimators. IEEE Trans. Biomed. Eng. 53, 43–52. doi: 10.1109/TBME.2005.859786

PubMed Abstract | CrossRef Full Text | Google Scholar

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, 82–87. doi: 10.1063/1.166141

PubMed Abstract | CrossRef Full Text | Google Scholar

Platisa, M. M., and Gal, V. (2006). Reflection of heart rate regulation on linear and nonlinear heart rate variability measures. Physiol. Meas. 27, 145–154. doi: 10.1088/0967-3334/27/2/005

PubMed Abstract | CrossRef Full Text | Google Scholar

Ponard, J. G., Kondratyev, A. A., and Kucera, J. P. (2007). Mechanisms of intrinsic beating variability in cardiac cell cultures and model pacemaker networks. Biophys. J. 92, 3734–3752. doi: 10.1529/biophysj.106.091892

PubMed Abstract | CrossRef Full Text | Google Scholar

Tan, C. O., Cohen, M. A., Eckberg, D. L., and Taylor, J. A. (2009). Fractal properties of human heart period variability: physiological and methodological implications. J. Physiol. 587(Pt 15), 3929–3941. doi: 10.1113/jphysiol.2009.169219

PubMed Abstract | CrossRef Full Text | Google Scholar

Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology. (1996). Heart rate variability. Standards of measurement, physiological interpretation, and clinical use. Eur. Heart J. 17, 354–381.

PubMed Abstract | Google Scholar

Tulppo, M. P., Hughson, R. L., Makikallio, T. H., Airaksinen, K. E., Seppanen, T., and Huikuri, H. V. (2001). Effects of exercise and passive head-up tilt on fractal and complexity properties of heart rate dynamics. Am. J. Physiol. Heart Circ. Physiol. 280, H1081–H1087.

PubMed Abstract | Google Scholar

Tulppo, M. P., Kiviniemi, A. M., Hautala, A. J., Kallio, M., Seppanen, T., Makikallio, T. H., et al. (2005). Physiological background of the loss of fractal heart rate dynamics. Circulation 112, 314–319. doi: 10.1161/CIRCULATIONAHA.104.523712

PubMed Abstract | CrossRef Full Text | Google Scholar

Viswanathan, G. M., Peng, C. K., Stanley, H. E., and Goldberger, A. L. (1997). Deviations from uniform power law scaling in nonstationary time series. Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics 55, 845–849. doi: 10.1103/PhysRevE.55.845

PubMed Abstract | CrossRef Full Text | Google Scholar

Wecht, J. M., Weir, J. P., and Bauman, W. A. (2006). Blunted heart rate response to vagal withdrawal in persons with tetraplegia. Clin. Auton. Res. 16, 378–383. doi: 10.1007/s10286-006-0367-y

PubMed Abstract | CrossRef Full Text | Google Scholar

West, C. R., Bellantoni, A., and Krassioukov, A. V. (2013). Cardiovascular function in individuals with incomplete spinal cord injury: a systematic review. Top. Spinal Cord Inj. Rehabil. 19, 267–278. doi: 10.1310/sci1904-267

PubMed Abstract | CrossRef Full Text

Willson, K., Francis, D. P., Wensel, R., Coats, A. J., and Parker, K. H. (2002). Relationship between detrended fluctuation analysis and spectral analysis of heart-rate variability. Physiol. Meas. 23, 385–401. doi: 10.1088/0967-3334/23/2/314

PubMed Abstract | CrossRef Full Text | Google Scholar