Heart Rate Fragmentation: A New Approach to the Analysis of Cardiac Interbeat Interval Dynamics

Background: Short-term heart rate variability (HRV) is most commonly attributed to physiologic vagal tone modulation. However, with aging and cardiovascular disease, the emergence of high short-term HRV, consistent with the breakdown of the neuroautonomic-electrophysiologic control system, may confound traditional HRV analysis. An apparent dynamical signature of such anomalous short-term HRV is frequent changes in heart rate acceleration sign, defined here as heart rate fragmentation. Objective: The aims were to: (1) introduce a set of metrics designed to probe the degree of sinus rhythm fragmentation; (2) test the hypothesis that the degree of fragmentation of heartbeat time series increases with the participants' age in a group of healthy subjects; (3) test the hypothesis that the heartbeat time series from patients with advanced coronary artery disease (CAD) are more fragmented than those from healthy subjects; and (4) compare the performance of the new fragmentation metrics with standard time and frequency domain measures of short-term HRV. Methods: We analyzed annotated, open-access Holter recordings (University of Rochester Holter Warehouse) from healthy subjects and patients with CAD using these newly introduced metrics of heart rate fragmentation, as well as standard time and frequency domain indices of short-term HRV, detrended fluctuation analysis and sample entropy. Results: The degree of fragmentation of cardiac interbeat interval time series increased significantly as a function of age in the healthy population as well as in patients with CAD. Fragmentation was higher for the patients with CAD than the healthy subjects. Heart rate fragmentation metrics outperformed traditional short-term HRV indices, as well as two widely used nonlinear measures, sample entropy and detrended fluctuation analysis short-term exponent, in distinguishing healthy subjects and patients with CAD. The same level of discrimination was obtained from the analysis of normal-to-normal sinus (NN) and cardiac interbeat interval (RR) time series. Conclusion: The fragmentation framework and accompanying metrics introduced here constitute a new way of assessing short-term HRV under free-running conditions, one which appears to overcome salient limitations of traditional HRV analysis. Fragmentation of sinus rhythm cadence may provide new dynamical biomarkers for probing the integrity of the neuroautonomic-electrophysiologic network controlling the heartbeat in health and disease.


INTRODUCTION
Heart rate variability (HRV) in healthy subjects, particularly over short time scales, is primarily attributable to fluctuations in vagal tone. The most recognizable manifestation of this parasympathetic influence is the oscillatory RR interval pattern (Figure 1) termed respiratory sinus arrhythmia (RSA) that results from the coupling between breathing and heart rate (Angelone and Coulter, 1964;Hirsch and Bishop, 1981;HRV, 1996;Stauss, 2003). However, beat-to-beat changes in the heart rate of healthy subjects not synchronized with respiration are also vagally mediated (Angelone and Coulter, 1964;Hirsch and Bishop, 1981). Therefore, a central interpretative framework underlying contemporary HRV analyses is one in which the degree of short-term variability of normal-to-normal (NN) sinus beats is used as a dynamical biomarker of cardiac vagal tone modulation (HRV, 1996;Billman, 2011).
This topic is of particular importance because parasympathetic regulation of sinus rhythm decreases with aging and organic heart disease (HRV, 1996;Kuo et al., 1999;Thayer et al., 2010). However, paradoxically, for some subjects in these high risk groups the amount of short-term variability actually increases (Figure 1). An apparent difference in the time series of vagally and non-vagally mediated HRV dynamics is their degree of smoothness, or conversely, their degree of fragmentation. Vagal tone modulation changes the heart rate in a progressive way. For example, with RSA, heart rate gradually increases and decreases with inspiration and expiration, respectively. When the coupling between heart rate and respiration is not as apparent, but the changes in heart rate are still driven by vagal tone modulation, the changes in heart rate are also gradual. In contrast, non-vagally mediated, shortterm heart rate variability has a distinct dynamical signature, namely more frequent changes in heart rate acceleration sign (Figure 1). In the "extreme" case of sinus alternans, the sign of heart rate acceleration changes every beat (Lewis, 1920;Geiger and Goerner, 1945;Friedman, 1956;Binkley et al., 1995). The presence of these abnormal variants of sinus rhythm limits the utility of traditional HRV analysis, since an increase in the overall amount of short-term variability can no longer be solely attributed to enhanced vagal tone modulation. Stein et al. (Domitrovich and Stein, 2002;Stein, 2002), coined the term "erratic sinus rhythm" to refer to prominent but apparently random variations in sinus cadence not attributable to vagal tone modulation and proposed a semi-quantitative approach to help identify them (Stein et al., 2005(Stein et al., , 2008. However, despite their association with increased cardiovascular risk and sick sinus syndrome (Bergfeldt and Haga, 2003), erratic sinus rhythm, sinus alternans, and their variants, have received scant clinical attention and the underlying mechanisms remain obscure.
As shown in Figure 1, the distinctions between the different classes of sinus arrhythmia may be difficult or impossible to discern from standard ECG recordings. The graphs of the NN interval time series and other representations of the data, such as Poincaré plots and Fourier spectra (not shown) may reveal clear differences in the structure of the fluctuations between physiologic and anomalous variability. However, in many cases, the differences are difficult to identify and especially to quantify.
These considerations led to the development of a novel approach to the analysis of short-term heart rate variability, termed heart rate fragmentation, accompanied by a set of simple-to-implement statistical metrics. A framework for the proposed approach is the concept that adaptive control of the heartbeat, particularly on short time scales, requires a hierarchy of interacting networks comprising neuroautonomic (especially the parasympathetic) and electrophysiologic components (sinus node pacemaker cells and their connections to the atrial syncytium). The integrity of these networks allows for their correlated function, evinced in part by the smoothness (fluency) of the output. At the same time, their functionality provides for sufficiently rapid (short-term or high frequency) responsiveness to physiologic stresses, while protecting against excessive volatility on a beat-to-beat basis.
A corollary concept is that network dysfunction, in general, and of the heart rate control system, in particular, is more likely to occur as the components of the network and their physiologic coupling start to break down. This degradative process should lead to increasing degrees of fragmentation. A key aspect of the fragmentation paradigm is that dysfunction or actual breakdown of one or more system components allows for the emergence of high frequency fluctuations that compete with or even exceed the shortest-term modulatory responsiveness of the vagal system. Therefore, a marker of this fragmentation on the surface ECG should be abrupt changes in the sign of heart rate acceleration, which may be periodic (as with classic sinus node alternans) or more random appearing (as with what has been termed "erratic sinus rhythm"). Such markers of fragmentation may be useful as correlates of cardiovascular aging and/or underlying organic heart disease.
Accordingly, we developed a set of fragmentation indices (see Methods) and applied them to beat-annotated, wellcharacterized 24-h Holter monitor recordings obtained from two very distinct clinical groups: healthy subjects and those with coronary artery disease (CAD). We analyzed three different time periods: the full day, putative awake and sleep periods. The primary hypotheses were that: (1) heart rate fragmentation would be higher in healthy old subjects than in younger ones for all FIGURE 1 | Examples of respiratory sinus arrhythmia and anomalous sinus rhythm. Electrocardiograms (Holter lead) from a healthy subject (top) and a patient with coronary artery disease (CAD) (middle), both from the present study. Normal-to-normal (NN) sinus interval time series from the healthy subject (bottom left) and of the patient with CAD (bottom right). The fluctuation patterns of the former time series are characteristic of phasic (respiratory) sinus arrhythmia, while that of the latter are indicative of an abnormal non-phasic sinus arrhythmia. To assist in visual comparisons, pale gray backgrounds are used for data from the healthy subject and light red for data from the patient with CAD, respectively. ECG voltage is given in arbitrary units (a.u). three time periods; and (2) heartbeat time series from patients with CAD would be more fragmented than those from healthy subjects. We also tested whether the fragmentation indices would outperform standard time and frequency domain measures, as well as nonlinear measures of short-term HRV in classifying heart rate time series from healthy subjects vs. those from patients with CAD.

Databases
We employed two long-term (∼24-h) ECG ambulatory databases from the Intercity Digital Electrocardiogram Alliance (IDEAL) study. The recordings are made available via the University of Rochester Telemetric and Holter ECG Warehouse (THEW) archives (http://thew-project.org/databases.htm). Putative waking and sleeping periods were estimated as the six consecutive hours of highest and lowest heart rates, respectively. These periods were calculated from the NN interval time series using a 6-h moving average window, shifted 15 min at a time.

HRV Analysis: Heart Rate Fragmentation Indices
From the ECG of each subject, the time series of the NN intervals, where t N i represents the time of occurrence of the i th normal sinus beat, and the time series of the differences between consecutive NN intervals (increments), The following four fragmentation indices were then computed from these time series: • The percentage of zero-crossing points in the increment time series, or equivalently, the percentage of inflection points (PIP) in the NN interval time series. (A t N i represents an inflection point if NN i × NN i + 1 ≤ 0, that is, if t N i is an instant of inversion of heart rate acceleration sign or of change to or from zero.) • The inverse of the average length of the acceleration/deceleration segments (IALS). An acceleration, deceleration segment is a sequence of NN intervals between consecutive inflection points for which the difference between two NN intervals is < 0 and > 0, respectively. The length of a segment is the number of NN intervals in that segment.
• The complement of the percentage of NN intervals in acceleration and deceleration segments with three or more NN intervals. This quantity is termed the percentage of short segments (PSS). • The percentage of NN intervals in alternation segments. An alternation segment is a sequence of at least four NN intervals, for which heart rate acceleration changes sign every beat. Such sequences follow an "ABAB" pattern, where "A" and "B" represent increments of opposite sign. This quantity is abbreviated, PAS.
By definition, the more fragmented a time series is, the higher the PIP, IALS, PSS, and PAS indices will be. We note that PAS quantifies the amount of a particular sub-type of fragmentation (alternation). A time series may be highly fragmented and have a small amount of alternation. However, all time series with large amount of alternation are highly fragmented. Given that the presence of non-sinus beats will increase fragmentation, we excluded segments encompassing non-sinus beats that started and ended at the inflection points preceding and following these non-sinus beats, respectively.
Finally, to assess the importance of beat annotation on fragmentation analyses, we also examined the full RR time series that include normal sinus beats as well as any supraventricular and ventricular ectopic beats.

HRV Analysis: Standard Measures
Standard techniques of HRV analysis are grouped into time and frequency (spectral) domain methods (HRV, 1996). A subset of the former, intended to quantify short-term variability, is based on the difference between consecutive normal-to-normal intervals ( NN, also termed NN increments); the latter on the spectral power of the NN intervals.
The following four traditional time and frequency HRV measures of short-term fluctuations were computed using opensource software (Goldberger et al., 2000) available at the PhysioNet website (www.physionet.org): • Time domain -pNNx measures: the percentage of NN > x ms. Here, we used x = 20 and 50 ms (Mietus et al., 2002). -rMSSD ("root mean square of successive differences"): square root of the mean of the squares of NN intervals. -SDSD ("standard deviation of successive differences"): standard deviation of the NN time series.
• Frequency domain -HF ("high frequency"): spectral power of the NN interval time series between 0.15 and 0.4 Hz.
These sets of time and frequency domain measures are widely interpreted to represent cardiac vagal tone modulation (HRV, 1996;Billman, 2011). By comparison, longer time scale fluctuations, are attributable to both sympathetic and parasympathetic influences (HRV, 1996;Thayer et al., 2010;Billman, 2013) and were, therefore, not considered here.

HRV Analysis: Non-linear Dynamical Indices
The following two widely used nonlinear short-term dynamical indices were computed: • Short-term detrended fluctuation analysis (DFA) exponent, α 1 . This measure (Peng et al., 1995) quantifies the correlations properties of a time series. The method is based on the assessment of the slope of the linear regression line of the log-log graph of F(n) vs. n. The function F(n) is the rootmean-square fluctuation of the integrated and detrended data, computed using windows of length n. For the analysis of heart rate time series, two indices, α 1 and α 2 , quantifying short and long-term behavior, respectively, have been proposed. Here, we focus on α 1 , which encompasses scales ranging from 4 to 11 beats, inclusively (Pikkujamsa et al., 1999). The correlation properties of time series with α ≃ 1.5 are similar to those of Brownian noise. In contrast, time series with α < 0.5 are anti-correlated. The former are smoother than the latter. • Sample entropy (SampEn). This measure (Richman and Moorman, 2000) quantifies the degree of irregularity of a signal. A higher SampEn value implies a more irregular, less predictable signal. Sample entropy is the negative of the natural logarithm of the conditional probability that the (m + 1) th components of two distinct segments match ( x i + m − x j + m < r) within the tolerance r, given that the first m components match within the same tolerance ( x i + l − x j + l < r, for 0 ≤ l ≤ m − 1).

Statistical Analysis
Spearman's rank and Pearson's product-moment correlation coefficients were used to quantify the dependence of: (i) the four novel indices of heart rate fragmentation, (ii) the traditional measures of short-term HRV, and iii) the two non-linear dynamical indices, short-term DFA exponent α 1 and SampEn, with the participants' age, using the THEW Healthy Subject Database. Statistical significance was set at a p-value <0.05. Logistic regression analysis methods were used to assess the relationships between presence of CAD and traditional, nonlinear and fragmentation indices in unadjusted models and models adjusted for age and gender. To facilitate comparisons among various HRV measures, we report normalized odds ratios (i.e., the odds ratio for a one standard deviation change in the measure).
The area under the receiver operating characteristic (AUC) curve was used to assess the goodness of fit of each model. The likelihood-ratio test was used to compare the goodness of fit of two nested models. All analyses were performed using raw measures except in the case of skewed variables whose logarithmic or quadratic transformation improved the models' goodness of fit. This improvement was only noted in the case of 24-h and daytime HF, 24-h and daytime SDSD and nighttime α 1 .

Changes in Heart Rate Dynamics with the Participants' Age in the Healthy Population
All four fragmentation indices significantly increased with the participants' age, for all three time periods, using either NN or RR interval time series (Table 1, Figure 2).
Traditional short-term HRV indices significantly decreased with the participants' age, for all three time periods.
The fractal α 1 exponent significantly increased with the participants' age during putative sleep time. For the other time periods, linear correlation analysis indicated an inverse relationship. However, in these cases, the Spearman coefficients were not significant.
Sample entropy significantly decreased with the participants' age during the putative wake and sleep periods. However, analyses of the 24-h period did not reveal any significant association between the two variables.
The percentage of supraventricular and ventricular premature beats significantly increased with the participants' age (Spearman r s = 0.27, p = 0.004 and r s = 0.29, p = 0.002, respectively).

Changes in Heart Rate Dynamics with Coronary Artery Disease
The values (median, 25 and 75th percentiles) of the new fragmentation indices for the groups of healthy subjects and patients with CAD, as well as of the traditional HRV and nonlinear indices, are presented in Table 2.
All fragmentation indices significantly (p < 0.0001) increased with the participants' age for all time periods in the group of patients with CAD, regardless of using NN or RR time series (Figure 2). The Pearson correlation coefficients varied between 0.250 and 0.529 for the NN time series and between 0.246 and 0.531 for the RR time series. The correlations for the 24-h and FIGURE 2 | Scatter plots of the traditional heart rate variability (rMSSD, pNN50, and HF), nonlinear (α 1 and SampEn) and fragmentation (PIP, IALS, PPS, and PAS) indices vs. the participants' age for the group of healthy subjects and patients with coronary artery disease (CAD), derived from the analysis of the full (∼ 24-h) period. The solid lines are the linear regression lines. rMSSD, root mean square of the successive differences; pNN50, percentage of differences between successive NN intervals above 50 ms; HF, high frequency spectral power; α 1 , detrended fluctuation analysis short-term exponent; SampEn, sample entropy; PIP, percentage of inflection points; IALS, inverse of the average length of the acceleration/deceleration segments; PPS, percentage of NN intervals in short segments; PAS, percentage of NN intervals in alternation segments. putative awake periods were relatively stronger than those for the sleep period.
Out of the 15 relationships tested, between each of the five traditional HRV measures and the participants' age, for each of the three time periods, only two were statistically significant. We found that pNN20 and pNN50 significantly decreased with the participants' age during the putative sleep period. Of the nonlinear indices, only DFA α 1 showed a significant association with participants' age. In this group, α 1 significantly decreased with the participants' age for all time periods.
A 1-year increase in age was associated with an increase of 14% in the odds of having CAD (odds ratio = 1.14, 95% confidence interval: 1.11-1.17, p < 0.0001). The AUC for the model with age as the only covariate was 0.853. Male sex carried a 3.54 fold increase in the odds of CAD (odds ratio = 3.54, 95% confidence interval: 2.17-5.78, p < 0.0001). The AUC for the null model with age and gender as the sole independent variables was 0.882.

Unadjusted analyses
In unadjusted analyses (Table 3), higher fragmentation indices were significantly associated with presence of CAD, for all time periods, using both NN and RR interval time series. Depending of the specific index and time period considered, a one-standard deviation increase in any of the fragmentation indices was associated with a 2.84-7.34-fold increase in the odds of CAD.
In comparison, the traditional short-term time and frequency domain HRV measures were inversely associated with presence of CAD for all time periods. However, only a subset of these measures, rMSSD and SDSD during sleep time, pNN50 during sleep and the 24-h period and pNN20 for all time periods, were significantly associated with CAD in unadjusted models. Of note, these models consistently performed worse than those with the fragmentation indices. For example, for pNN20, the best performing of the HRV measures, a one-standard deviation increase in the value of this variable was only associated with a 26, 77, and 44% increase in the odds of CAD, for the awake, sleep and 24-h periods, respectively.
Lower values of α 1 , for the awake and 24-h periods, and of SampEn for the sleep period were also significantly associated with presence of CAD. Overall, α 1 was a stronger correlate of CAD than traditional HRV measures, but not as strong as the fragmentation indices. SampEn, even for the sleep period, was among the weakest correlates of CAD. Figure 3 shows the normalized histograms of the traditional heart rate variability (rMSSD, pNN50, and HF), nonlinear (α 1 and SampEn) and fragmentation (PIP, IALS, PPS, and PAS) indices for the groups of healthy subjects and patients with coronary artery disease, for the 24-h period.

Adjusted Analyses
Fragmentation indices remained positively associated with CAD in models adjusted for age and sex (Table 4). Furthermore, the models with any of these indices fitted the data better than the ones with only age and sex, for all time periods, regardless of whether NN or RR time series were used. FIGURE 3 | Normalized histograms of the traditional heart rate variability (rMSSD, pNN50, and HF), nonlinear (α 1 and SampEn) and fragmentation (PIP, IALS, PPS, and PAS) indices for the groups of healthy subjects (blue) and patients with coronary artery disease (red), for the 24-h period. rMSSD, root mean square of the successive differences; pNN50, percentage of differences between successive NN intervals above 50 ms; HF, high frequency spectral power; α 1 , detrended fluctuation analysis short-term exponent; SampEn, sample entropy; PIP, percentage of inflection points; IALS, inverse of the average length of the acceleration/deceleration segments; PPS, percentage of NN intervals in short segments; PAS, percentage of NN intervals in alternation segments. Adding any of the fragmentation indices derived from NN interval time series to a model of CAD with the percentages of supraventricular premature beats (% SVPBs) significantly increased its performance (p < 0.00001) for all time periods. Specifically, while the AUC for the model with the % SVPBs was 0.754, the AUCs for the models that also included PIP, IALS, PPS, or PAS derived from NN intervals time series, were 0. 852, 0.856, 0.866, and 0.771, respectively. None of the associations between traditional time and frequency domain measures and CAD remained significant for any of the time periods when the models were adjusted for age and sex (Table 4). For the nonlinear measure α 1 , the associations remained significant during the 24-h and awake periods. For SampEn, the association with CAD during awake was significant, while the relationship during sleep time lost significance. The adjusted models with α 1 , for the 24-h and awake periods, and with SampEn, for the awake period, were superior to the model with only age and sex.

Relationship between HRV Measures and Mean Heart Rate
In both healthy subjects and those with CAD, the fragmentation indices, PIP, IASL, and PSS, were not significantly correlated with 24-h mean heart rate (Figure 4). The same was true for the group of patients with CAD. PAS was weakly correlated with mean heart rate in both healthy subjects and patients with CAD. The correlation was positive, r = 0.221 (p = 0.021), for the group of healthy subjects and negative, r = −0.146 (p = 0.020) in the group of those with CAD.
In contrast, the traditional time and frequency domain measures were strongly correlated with mean heart rate. The standardized Pearson correlation coefficients varied between 0.501 and 0.676 (p < 0.0001) in the group of healthy subjects and between 0.111 and 0.574 in the group of patients with CAD (p < 0.0001 for all correlations, but the one with SDSD).
SampEn was only weakly correlated with mean heart rate. The correlation was negative in the healthy group, r = −0.185, p = 0.054 and positive in the CAD group, r = 0.287, p < 0.0001. α 1 was positively correlated with mean heart rate in patients with CAD r = 0.156, p < 0.013. For the healthy subjects, the correlation was not significant.

DISCUSSION
This study is of potential interest because it presents a new way of assessing short-term HRV under free-running (spontaneous) conditions. The novel methodology and findings are described under the rubric of sinus rate fragmentation (or conversely, smoothness or "fluency"). The conceptual framework is described in the Introduction. We found that the degree of fragmentation of the NN and RR time series, derived from 24-h Holter monitoring, varied directly as a function of cross-sectional age in a cohort of healthy young to elderly male and female subjects in sinus rhythm. In this group, older age was associated with increased fragmentation. This correlation was noted for the entire 24-h period, as well as, during putative wake and sleep periods. Furthermore, we found that the fragmentation indices outperformed standard time and frequency domain measures, as well as, two widely used nonlinear measures (DFA α 1 and SampEn), in separating healthy subjects from patients with CAD.
FIGURE 4 | Scatter plots of the traditional heart rate variability (rMSSD, pNN50, and HF), nonlinear (α 1 and SampEn) and fragmentation (PIP, IALS, PPS, and PAS) indices vs. mean heart rate (in beats per minute, bpm) for the group of healthy subjects (blue dots) and those with coronary artery disease (CAD, red circles), derived from the analysis of 24-h NN interval time series. The blue and red lines are the linear regression lines for the healthy and CAD groups, respectively. rMSSD, root mean square of the successive differences; pNN50, percentage of differences between successive NN intervals above 50 ms; HF, high frequency spectral power; α 1 , detrended fluctuation analysis short-term exponent; SampEn, sample entropy; PIP, percentage of inflection points; IALS, inverse of the average length of the acceleration/deceleration segments; PPS, percentage of NN intervals in short segments; PAS, percentage of NN intervals in alternation segments.
In this study, we chose to analyze open access Holter data from groups of subjects whose clinical status was well-characterized and presented very sharp population differences: a group of healthy subjects and a group of patients with overt CAD. The healthy subjects were, on average, 20 years younger than the patients with CAD. Therefore, these two group were robustly separated by a model that simply incorporated age and gender (AUC = about 0.88). We expected the combination of older age and overt cardiovascular disease in the CAD group to enhance the ability of quantitative methods of HRV to unambiguously discriminate between patients and healthy individuals.
A potential link between aging and a variety of overt cardiovascular disease processes is the role of inflammation and fibrosis (Biernacka and Frangogiannis, 2011;Ghiassian et al., 2016), which decrease the amount and/or effectiveness of physiologic vagal tone modulation and promote the breakdown of regulatory networks, such as those controlling heart rate. Therefore, all short-term HRV measures were expected to change in the same directional way with aging and cardiovascular disease, for all time periods. In this regard, and in accord with "canonical" HRV precepts (HRV, 1996;Stauss, 2003;Billman, 2011), we hypothesized that traditional short-term HRV measures would decrease with cross-sectional age in the healthy group and that these measures would be lower for the patients with CAD than healthy subjects. In addition, we hypothesized that the fragmentation indices would increase with the participants' age in the healthy group and that they would be higher for patients with CAD than healthy subjects.
We found strong correlations between traditional HRV measures and the participants' age in the healthy group, for all time periods. These findings were in agreement with the generally accepted idea that HRV measure are useful to assess changes in heart rate dynamics with healthy aging (Pikkujamsa et al., 1999). In contrast, for the group of patients with CAD, traditional HRV indices, with only two exceptions, pNN20 and pNN50 during putative sleep, did not significantly change with the participants' age. This counterintuitive finding (Figure 2) may be due to the confounding effects of heart rate fragmentation, which may increase high-frequency variability not due to physiologic respiratory-vagal modulation (Domitrovich and Stein, 2002;Stein, 2002;Stein et al., 2005Stein et al., , 2008. Furthermore, the ability of traditional short-term HRV measures to separate the healthy and CAD groups was also surprisingly poor. In models adjusted for age and sex, none of the traditional HRV measures significantly discriminated these two groups. Even in unadjusted models, the discriminatory power of conventional HRV measures was not consistent across time periods. In particular, HF power, traditionally interpreted as a measure of vagal tone modulation, did not discriminate the two groups for any of the time periods considered. Our results are consistent with previous cautionary reports about the utility of traditional HRV measures to assess vagal tone modulation especially with advanced age or underlying heart disease (Stein, 2002;Stein et al., 2005;Burr, 2007;Billman, 2013).
The nonlinear indices also did not provide consistent results. For example, α 1 significantly increased with the participants' age during sleep, in both Pearson and Spearman correlation analyses.
However, an inverse relationship was found for the awake and 24h periods, in Pearson but not in Spearman analyses. In addition, lower α 1 values were significantly associated the presence of CAD during the awake and 24-h periods, but not during sleep. The degree of randomness of heart rate time series, measured by SampEn, significantly decreased with the participants' age during the awake and sleep periods. Furthermore, while in an unadjusted model, a one standard deviation increase in the degree of randomness of heart rate time series was associated with a 36% decrease in the odds of CAD during sleep time, in a models adjusted for age and sex, a one standard deviation increase in SampEn was associated with a 57% increase in the odds of CAD during awake period. For the other time periods, the associations were not significant.
The inconsistencies of the two nonlinear methods, α 1 and SampEn, are not entirely unexpected. For example, the fluctuation function, F(n), in DFA analysis of 24-h NN interval time series usually presents a crossover separating "short-term" from "longer-term" behavior. However, the scale at which that crossover occurs may vary substantially from subject to subject. In addition, the degree of linearity of F(n) also tends to vary from subject to subject. SampEn, on the other hand, can be affected by nonstationarities that are common in real world data. In addition, fragmented time series can be highly predictable (leading to low SampEn), as in the case of those with a high density of alternation, or highly irregular, as in the case of "erratic" sinus rhythm (leading to increased SampEn).
The limitations of traditional and newer HRV methods, exemplified by the results reported above and those described by other investigators (Stein, 2002;Stein et al., 2005;Burr, 2007;Billman, 2013), help motivate the on-going searches for alternative approaches. The introduction of the concept of fragmentation of heart rate dynamics, accompanied by a set of metrics for its quantification, are part of this exploration.
Speculatively, possible mechanisms of the observed fragmentation include the breakdown of one or more components of the regulatory network controlling heart rate dynamics. An obvious first question would be whether the higher fragmentation values in CAD vs. the healthy group could simply be due to SVPBs mislabeled as normal sinus beats. While the THEW website describes that three lead Holter monitor recordings were first processed using an automated beat annotation program and then subjected to visual review and adjudication, the possibility that some of the beats labeled as N are actually subtle SVPBs, and not sinus beats, cannot be absolutely excluded. To address this possible confounder, one would assume that the recordings with the highest likelihood of containing hidden SVPBs would be those with the highest percentage of labeled SVPBs. Therefore, to assess whether our results were likely a consequence of mislabeled SVPBs, we compared the performance of a model of CAD with the % SVPBs, as the sole independent variable to that of a model with the % SVPBs and each fragmentation index. We found that all fragmentation indices added significant information to the model with % SVPBs. This finding supports the contention that fragmentation is not a surrogate measure of "hidden" supraventricular ectopy.
A second question would be whether abnormalities in breathing dynamics could be responsible for the fragmentation of heart rate dynamics through respiratory-cardiac coupling. Since we did not have a direct measure of respiration, we cannot exclude the possibility that inter-breath interval time series were themselves fragmented. However, such a mechanism is unlikely since beat-to-beat changes in the sign of heart rate acceleration are above the frequency response of the vagalsinus modulatory system. In fact, the coupling between the sinus node and the vagus tends to drop-off at high respiratory rates (Angelone and Coulter, 1964;Hirsch and Bishop, 1981). Furthermore, the most erratic variants of sinus rhythm are seen in the populations of older healthy subjects and those with organic heart disease (Stein et al., 2005(Stein et al., , 2008, groups with the most impaired vagal modulatory capacity and, therefore, those least likely to show very high frequency coupling between autonomic and electrophysiologic components. The specific electrophysiologic bases for fragmentation of heart rate dynamics remain to be determined. More than one mechanism may be contributory. For example, alternans phenotypes could be due to sinus node exit block or to very subtle atrial bigeminy with SVPBs originating near or even within the sino-atrial (SA) node (Geiger and Goerner, 1945). Another mechanism that could account for fragmentation would be modulated sinus node parasystole (Jalife et al., 1986), an arrhythmia in which two pacemaker sites in the SA area show bidirectional coupling and appear to "compete" for control of the heartbeat. Under certain parameter regimes, such coupling may induce a variety of alternating NN patterns (Jalife et al., 1986). The underlying electrophysiologic mechanisms to account for fragmentation may also involve perturbations of internal pacemaker "clocks" in the SA node (Lakatta et al., 2010). From a pathophysiologic viewpoint, mechanisms related to altered conduction and/or abnormal automaticity all reflect instabilities in the parasympathetic-SA node-atrial network. Given this substrate of instability, highly fragmented (whether erratic or periodic) types of NN patterns may represent pre-or even proarrhythmic markers. Our findings that fragmentation is increased in the elderly and in those with established CAD support this contention. Fragmentation would be of high interest if it were a forerunner of arrhythmias such as atrial fibrillation or other tachyarrhythmias in which the control network becomes so unstable that sinus node function is overridden by ectopic stimuli. Whether fragmentation is an independent risk marker of cardiovascular mortality related to heart failure also remains to be determined, as does any relationship to classical sick sinus syndrome.
More generally, the findings here support a modification in the standard classification of sinus rhythm into "phasic" and "non-phasic" variants (Faulkner, 1930;Hirsch and Bishop, 1981;Fisch and Knoebel, 2000). The first category refers to the oscillations in heart rate that are coherent with respiration and are most marked in younger individuals at rest, during deep sleep or with meditation (classic RSA). Non-phasic sinus arrhythmia, a term that has largely disappeared from the clinical lexicon, has been used to refer to a variety of sinus variants without this strict periodicity, including erratic sinus rhythm, and usually connotes abnormal sinus function (Stein et al., 2005). However, non-phasic types of sinus arrhythmia may also occur as physiologic variants, e.g., during exercise and recovery.
An alternative schema would be classify sinus rhythm into phasic and non-phasic types, and then sub-divide non-phasic into either physiologic due to short term trends but without tight respiratory coupling and non-physiologic, i.e., fragmented categories. However, we emphasize that fragmentation analysis per se does not separate phasic and non-phasic variants into two discrete bins. Rather, it quantifies, in a continuous way, the degree to which fragmentation is present.
Heart rate fragmentation may account for some of the abnormal patterns in Poincaré and other maps previously reported (Woo et al., 1992;Brouwer et al., 1996;Huikuri et al., 1996;Domitrovich and Stein, 2002;Stein et al., 2005Stein et al., , 2008Gladuli et al., 2011;Makowiec et al., 2015). Such maps contain important information about the temporal structure of a time series. However, they are difficult to quantify in a physiologically interpretable way. Commonly employed metrics such as SD1, SD2, and SD1/SD2, only measure linear properties of the data that are also captured by time domain HRV measures such as rMSSD and SDSD (Brennan et al., 2001). If heart rate fragmentation is found to be one of the mechanisms underlying such abnormal patterns, the metrics introduced here may help identify, in a fully automated way, the time series associated with certain types of anomalous Poincaré plots.
Other fragmentation-related indices may also prove useful. Examples include the densities of: (i) "hard edges, " defined as inflection points for which NN i × NN i + 1 < 0; (ii) "soft edges, " defined as NN i × NN i + 1 = 0, where NN i = NN i + 1 ; (iii) "short segments, " defined as acceleration/deceleration segments encapsulated between "hard edges" or "soft edges;" and (iv) segments for which heart rate does not change. Additionally, symbolic dynamical analysis of heart rate increment time series, where words are analyzed in terms of the number of edges they contain may also prove useful in this context. Finally, the fragmentation indices have a number of attractive features. First, these indices are independent of the mean heart rate (Figure 4). The only exception is PAS, which is not a general fragmentation index, but quantifies a particular type of fragmentation (pattern of the type "ABAB, " where "A" and "B" represent increments of opposite sign). In contrast, traditional short-term time and frequency domain measures showed highly significant negative associations with mean heart rate, both in the group of healthy subjects and of those with CAD. These results are in line with those reported in other studies (Monfredi et al., 2014;Sacha, 2014). Second, by construction, the fragmentation indices (including PAS) are also independent of the amplitude of the time series. This feature is due to the fact that accelerations/decelerations were defined as increments/decrements in heart rate of any magnitude. Thus, two time series with fluctuation patterns that only differ in amplitude (e.g., time series, u i and v i , for which u i /v i = c, where c is a constant) will have exactly the same degree of fragmentation. Future studies will be needed to explore whether the use of a threshold > 0 in the definition of accelerations and decelerations further increases the discriminatory power of these measures in HRV analyses. Third, the fragmentation indices are among the measures least affected by nonstationarities. The reason is that the operation of calculating the increment time series, used to detect the inflection points (i.e., changes in heart rate acceleration sign), detrends the data. Fourth, the fragmentation indices can be computed using NN or RR interval time series. Indeed the use of the latter did not impair the discriminatory power of the fragmentation indices for the populations studied here. This finding, if validated, may facilitate fully automated implementations of the method. Future studies will also help determine the effect of data length on the confidence intervals of the fragmentation indices.

CONCLUSION
Analysis of short-term HRV is enhanced by a set of computational tools that quantify the fragmentation of heartbeat variability, defined by abrupt changes in the sign of HR acceleration. In a Holter monitor database from healthy subjects, the degree of fragmentation increased with the participants' age.
Furthermore, fragmentation measures outperformed traditional short-term measures of HRV in discriminating a group of patients with CAD and from the healthy subjects. Fragmentation of sinus rhythm cadence may support a new class of dynamical biomarkers that probe the integrity of the regulatory network comprising neuroautonomic, sinus node and atrial components.

AUTHOR CONTRIBUTIONS
MC and AG developed the fragmentation concept and method. RD directed the statistical analysis. All three authors contributed to the interpretation of the findings and worked collaboratively on the manuscript.