Heart Rate Fragmentation as a Novel Biomarker of Adverse Cardiovascular Events: The Multi-Ethnic Study of Atherosclerosis

Background: A major objective of precision medicine is the elucidation of non-invasive biomarkers of cardiovascular (CV) risk. Recently, we introduced a new dynamical marker of sino-atrial instability, termed heart rate fragmentation (HRF), which outperformed traditional and nonlinear heart rate variability metrics in separating ostensibly healthy subjects from patients with coronary artery disease. Accordingly, we hypothesized that HRF may be a dynamical biomarker of adverse cardiovascular events (CVEs). Methods: This study employed data from a cohort of participants in the Multi-Ethnic Study of Atherosclerosis (MESA), a prospective study of sub-clinical heart disease. Interbeat interval time series (n = 1963), derived from the electrocardiographic channel of the polysomnogram study, were analyzed using the newly introduced metrics of fragmentation, as well as traditional heart rate variability (HRV) indices and the short-term detrended fluctuation analysis exponent. Cox regression analysis was used to assess the association between HR dynamic indices and CV outcomes in unadjusted and adjusted models. Results: The mean (± SD) follow-up time was 2.97 ± 0.63 years. In adjusted models, higher fragmentation was significantly associated with incident CVEs (number of events; hazard ratio [95% confidence interval]: n = 72, 1.43 [1.16–1.76]) and CV death (n = 21; 1.65 [1.15–2.36]). The traditional HRV and the fractal indices were not associated with CVEs or CV death. The most discriminatory fragmentation indices added significant value to Framingham and MESA CV risk indices in all analyses. Conclusion: Our findings show that HRF has promise as a non-invasive, automatable biomarker of CV risk. The basic mechanisms underlying fragmentation remain to be delineated. Its association with incident outcomes raises the possibility of connections to degenerative changes in the multisystem network controlling SAN function.


INTRODUCTION
This study describes a novel noninvasive biomarker of cardiovascular (CV) risk based on heart rate dynamics. In healthy adults at rest and during sleep, the highest frequency at which the sino-atrial node (SAN) rate fluctuates varies between ∼0.15 and 0.40 Hz (Figures 1A 1 -A 4 ). These oscillations, referred to as respiratory sinus arrhythmia, are due to vagally-mediated coupling between the SAN and breathing. However, not all fluctuations in heart rate (HR) at or above the respiratory frequency are attributable to vagal tone modulation. Under pathologic conditions, an increased density of reversals in HR acceleration sign, not consistent with short-term parasympathetic control, can be observed ( Figures 1B 1 -B 4 ). This dynamical biomarker of electrophysiologic instability has recently been identified and termed heart rate fragmentation (HRF) (Costa et al., 2017a). A set of metrics (computational probes) for its quantification was also introduced (Costa et al., 2017a,b).
Perhaps the most explicit example of HRF is the subtle supraventricular arrhythmia termed sinus node alternans (Binkley et al., 1995), in which the time between consecutive sinus beats oscillates between two values, short (S) and long (L) following an SLSL pattern. However, HRF includes not only pure (2:1) sinus node alternans but also quasi-periodic and more irregular variants of normal-to normal (NN) alternation. As Figures 1A 3 ,B 3 illustrates, clinical recognition of such patterns is difficult from standard electrocardiograms (ECGs). Traditional heart rate variability (HRV) indices and other metrics such as those derived from Poincaré plots (Figures 1A 4 ,B 4 ,C) may be also of limited use in the identification of HRF patterns. The basic mechanisms of fragmentation, involving either anomalous sinus beats (Lewis, 1920) or supraventricular ones originating near the SAN, are unresolved (Costa et al., 2017a,b).
The potential importance of HRF is several-fold. First, it produces a high degree of short-term variability that may be mistaken as a marker of healthy vagal control when standard measures of short-term HRV are used. Second, its presence supports the delineation of a new class of biomarkers of cardiac risk. The latter is premised on the conjectured link between HRF and the breakdown in one or more components of the control system (and/or in their interactions) regulating SAN function. Notably, earlier reports of what were first termed "sinus extrasystoles" as well as sinus alternans (Lewis, 1920;Binkley et al., 1995) were from patients who were older or had organic heart disease. Third, the investigation of fragmentation may yield new insights into SAN functionality in health, aging and disease.
In recent studies (Costa et al., 2017a,b) we analyzed annotated Holter recordings [University of Rochester Telemetric Holter ECG Warehouse (THEW)] from healthy subjects and patients with advanced coronary artery disease (CAD) using the newly devised HRF metrics. Fragmentation was found to significantly increase as a function of the participants' age in both the healthy population and the one with CAD. In contrast, most short-term HRV indices did not significantly change with the participants' age in the CAD group. Furthermore, fragmentation was higher in patients with CAD than in healthy subjects, during both estimated awake and sleep periods, while traditional HRV metrics did not discriminate the two groups.
The general motivation for the present study was to assess the potential utility of the novel indices of HRF as predictors of adverse cardiovascular events (CVEs) and CV mortality, using the large Multi-Ethnic Study of Atherosclerosis (MESA). This ongoing prospective cohort study (Bild et al., 2002) was designed to investigate the prevalence, correlates and progression of subclinical cardiovascular disease (CVD) in a multi-ethnic population free of overt clinical CVD at study entry. We specifically hypothesized that HRF would: (1) be positively associated with cross-sectional age; (2) be positively associated with incident CVEs and CV death; and (3) outperform traditional HR dynamical measures. In addition, we also sought to determine if fragmentation metrics added value to prediction tools computed from clinical measures, namely the Framingham Heart Study (D'Agostino et al., 2008) and MESA CV risk indices (McClelland et al., 2015).

Study Population and Data Collection
The MESA study has been described in detail previously (Bild et al., 2002). Briefly, over a period of approximately 2-years, starting in July 2000, 6,814 persons between the ages of 45 and Frontiers in Physiology | www.frontiersin.org 84 years of age without evident clinical CVD were recruited at six field centers in the US. Institutional review boards from each study site approved the conduct of this study, and written informed consent was obtained from all participants. A sleep ancillary study was conducted in conjunction with MESA's fifth examination (2010)(2011)(2012)(2013). The study enrolled 2,060 participants who underwent unattended, inhome polysomnography (PSG) following a standardized protocol (Redline et al., 1998). The data obtained using the 15-channel Compumedics Somte System (Compumedics LTd., Abbottsville, Australia) were scored at the Brigham and Women's Hospital centralized reading center by trained technicians using published guidelines (Redline et al., 2007). The apnea-hypopnea index (AHI) was calculated based on the average number per hour of sleep of all apneas plus hypopneas associated with ≥ 3% oxygen desaturation or arousal.
The ECG channels, sampled at 256 Hz, were processed using Compumedics Somte software for detection and classification of the QRS complexes (R-points) as normal sinus, supraventricular premature or ventricular premature complexes. The automated annotations were reviewed by a trained technician, who made appropriate corrections. Both the NN and the R-to-R (RR) interval time series were analyzed in the present study.
Participants with one or more of the following were excluded: poor signal quality (n = 35), pacemaker (n = 13), atrial fibrillation (AF) at the time of the PSG (n = 22), < 2 h of combined sleep periods scored as rapid eye movement (REM), stage 1, 2, 3, or 4 (n = 16), and < 75% normal sinus beats between sleep onset and termination (n = 11). Participants with CVEs before the PSG (n = 185) and seven others for whom the last recorded followup was prior to the PSG were excluded from the analyses of the associations between HR dynamical metrics and incident CVEs. These participants were included in analyses of CV mortality.

Clinical Follow-Up and Event Classification
In addition to clinical exams, participants are followed every 9-12 months to inquire about hospital admissions, CV outpatient diagnoses and procedures, and deaths. Discharge diagnosis codes are obtained for all hospitalizations and medical records are obtained when heart failure, myocardial infarction, stroke, or death are reported. For those over age 65 and enrolled in fee-for-service Medicare, claims data are also used to identify diagnosis and procedure codes. Trained personnel abstract any hospital records suggesting possible CVEs, which are then adjudicated by physicians at the coordinating center. Nonfatal endpoints in MESA include congestive heart failure, angina, myocardial infarction, percutaneous coronary intervention, coronary bypass grafting or other revascularization procedure, resuscitated cardiac arrest, peripheral arterial disease, stroke (non-hemorrhagic) and transient ischemia attack (TIA). Cardiovascular deaths, as adjudicated by committee review, included fatalities directly related to stroke or coronary heart disease. For other deaths, the underlying cause are obtained through state or vital statistics departments. The definition and adjudication of these events have been described in detail previously (Bild et al., 2002;Bluemke et al., 2008;Yeboah et al., 2012). The cut-off date for the surveillance period was December 31, 2014.

Fragmentation Analysis
Fragmentation analysis was performed for 1963 subjects using both NN and RR interval time series. Fragmentation analysis is described in detail in Costa et al. (2017a,b). Briefly, original interbeat interval time series, {s i }, 1 ≤ i ≤ L (L, time series length) were mapped to a ternary symbolic sequence as follows: "−1" if NN i < −4 ms, "0" if −4 ≤ NN i ≤ 4 ms, and "1" if NN i > 4 ms (Costa et al., 2017b). Note that, since the ECG signals were sampled at 256 Hz, the resolution of the interbeat interval time series is 1/256 ∼ 4 ms. Therefore, only NN (or RR) intervals whose difference was > 4 ms or < −4 ms were considered different from each other.
To assess the prevalence of dynamical patterns with increasing degrees of fragmentation, the percentages of sequences of 4 consecutive symbols, w i = (s i , s i+1 , ..., s i+l−1 ), 1 ≤ i ≤ L − l + 1, termed "words, " with 0, 1, 2, and 3 inflection points were calculated. We refer to these word classes as W 0 , W 1 , W 2 , and W 3 , respectively. The full lexicon that comprises 81 different words is given in Costa et al. (2017b). Words derived from the NN (RR) interval time series were termed NN (RR) words. The words in groups W 0 and W 1 are the least fragmented (most "fluent"), those in groups W 2 and W 3 are the most fragmented.
All fragmentation metrics were calculated for the entire sleep period. In the case of PIP, standard sleep stage periods (awake, rapid eye movement [REM], stages N1, N2, and N3), in addition to periods of awake before sleep onset or after sleep termination were also analyzed. Individuals were included if cumulatively they had at least 2000 NN intervals during a given sleep stage. Thus, the number of participants in each of these sub-analyses was a fraction of the number of participants in the analyses of the entire sleep period. Specifically, the number of participants (the number of those with incident CVEs) for the different periods was: awake, 1379 (61); REM, 1451 (54); N1, 1056 (51); N2, 1695 (70); N3, 791 (23); and awake before sleep onset or after sleep termination, 1449 (58).

Traditional HRV Analysis
The following traditional time domain HRV indices (HRV, 1996) were calculated for 1963 subjects using NN interval time series: (1) the average of all NN intervals (AVNN), (2) mean of the standard deviations (SDs) of NN intervals in all 5-min segments (SDNNIDX), (3) the square root of the mean of the squares of differences between adjacent NN intervals (rMSSD), and (4) the percentage of differences between adjacent NN intervals that are greater than 50 ms (pNN50). The following traditional frequency domain HRV indices were calculated: (1) the total spectral power of all NN intervals between 0.15 and 0.4 Hz (HF) and (2) the ratio of low to high frequency power (LF/HF). Each of these metrics was calculated using a 5-min sliding window (without overlap), with more than 150 beats and more than 75% NN intervals, between sleep onset and sleep termination. Power spectrum estimates were obtained using the Lomb periodogram method, which does not require missing data points (due to removal of ectopic, misdetections, and artifact) in a time series to be interpolated (Moody, 1992(Moody, , 1993. A total of 170,527 windows were analyzed. For each subject, the values from the different windows were averaged. The source code used for these computations is published on The Research Resource for Complex Signals, (PhysioNet) website (Goldberger et al., 2000;Mietus and Goldberger, 2008).

Detrended Fluctuation Analysis: Short-Term Scaling Index
Detrended fluctuation analysis (Peng et al., 1995) was developed to quantify the correlation properties of a time series. The methods is based on the assessment of the slope of the regression line of the logarithm of F(n) vs. the logarithm of n, where F(n) is the root-mean-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. We focused on α 1 , defined as the slope of F(n) vs. n, for 4 ≤ n ≤ 11 (Pikkujamsa et al., 1999). Time series for which α ∼ 0.5 are uncorrelated (random). Those with α > 0.5 and those with α < 0.5 are long-range correlated and anti-correlated, respectively.
Discontinuities in the NN interval time series due to premature beats, misdetections and artifact were dealt with in the following way. If the gap was < 3 s (typically due to the removal of an ectopic beat), an interpolated beat was inserted. If the gap was wider (typically due to ECG artifact that prevents accurate detection of the peak of the QRS complexes), then the segments that preceded and followed the gap, were "stitched" together. The α 1 exponent was calculated for non-overlapping segments of 1,000 intervals and the results were averaged. The source code used in these computations is published on the PhysioNet website (Goldberger et al., 2000;Mietus et al., 2001).

Statistical Analysis
Continuous variables are summarized as median, first and third quartiles, unless otherwise indicated. Categorical variables are presented as numbers and percentages.
The associations between independent variables and both incident CVEs and CV mortality were assessed using Cox proportional hazard analysis. Efron's method (Efron, 1977) was used to handle ties. Failure time in the individuals with incident CVEs was the time between the PSG study and the time of event diagnosis. For participants without CVEs, the failure time was the time between the PSG study and the latest follow-up, death, or loss to follow-up. Statistical significance was set at a pvalue <0.05. The independent variables were: the fragmentation indices, PIP, W 0 , W 1 , W 2 , and W 3 , derived from both NN and RR interval time series; the HR dynamical indices: AVNN, SDNNIDX, rMSSD, pNN50, HF, and HF/LF; and the short-term fractal index DFA α 1 .
Standardized hazard ratios (per one-SD increase in the independent variable) were calculated with associated 95% confidence intervals (CI). The assumption of proportional hazards was tested using a global test based on Schoenfeld residuals (Grambsch and Therneau, 1994). No violations were noted. The predictive power of the survival models was assessed using Harrell's C statistic. The likelihood ratio test was used to compare the fit of two nested models (null model vs. null model + HR dynamical metric). The three null models were models 2, 3, and 4 described above. The null hypothesis for each of the likelihood ratio tests was that the two nested models fitted the data equally well. Rejection of the null hypothesis implied that the larger model fitted the data better, indicating that a given HR dynamical metric added value to the null model.
Linear regression models with a quadratic term were used to describe the possible nonlinear (U-shaped) relationship between: i) HRF and short-term HRV indices (example of a model with PIP and rMSSD, PIP = β 1 * ln(rMSSD) + β 2 * [ln(rMSSD)] 2 + α, where α is a constant); and ii) short-term HRV indices and the participants' age.
In all analyses, the variables W 0 , rMSSD, pNN50, SDNN, or HF were logarithmically transformed since the models with these log-transformed variables fitted the data better than those with untransformed ones. In all other cases, the analyses were performed using untransformed variables.
Characteristics of MESA participants without and with a CVE during follow-up are summarized in Table 1. Individuals who developed CVEs were older and more likely to be male and have diabetes. They tended to have higher seated HR and higher systolic blood pressure. In addition, this risk group tended to have lower sleep efficiency and a higher apnea-hypopnea index. The differences between those who died and did not were qualitatively similar to the differences between those with and without incident CVEs.

Relationship of HRF, Traditional HRV and Nonlinear Indices With the Participants' Age
All fragmentation indices, derived from either NN or RR interval time series, were significantly associated with the participants' age. The Pearson correlation coefficients (ρ) for PIP, W 0 , W 1 , W 2 , and W 3 were 0.35, −0.17, −0.31, 0.10, and 0.35, respectively. PIP (Figure 2) and the percentages of fragmented words W 2 and W 3 increased with the participants' age at the estimated rates of 0.28, 0.09, and 0.35% per year of age, respectively. The percentages of fluent words W 0 and W 1 decreased with the participants' age at the rates of −0.06 and −0.39% per year, respectively. Slightly higher rates of change were observed for the indices derived from RR interval time series (not shown). The short-term traditional HRV indices, rMSSD, pNN50 and HF, did not vary linearly with the participants' age. Instead, the association between age and short-term variability depended on the participants' age itself. Figure 2 illustrates, using the representative example of rMSSD, how the amount of shortterm variability varied across different age groups. Variability was higher in both the lowest (<55 years) and highest (≥85 years) age groups compared to intermediate ones (U-shape relationship).
The short-term fractal index, α 1 showed a small but significant decrease with cross-sectional age. The Spearman correlation coefficient was −0.18 (p = 0.000).

Unadjusted Analyses of Risk of Incident CVEs
All fragmentation indices, calculated from both the NN and the RR time series, were significantly associated with incident events ( Table 2: Model 1). The association was positive for fragmented words, W 2 and W 3 , as well as PIP, and negative for fluent (less fragmented) words, W 0 and W 1 . The most discriminatory of the fragmentation indices was the word W 1 derived from the RR interval time series. A one-SD increase in W 1 was associated with a 44% (95% CI: 28-66%) decrease in the rate of CVEs. This variable performed comparably to the Framingham Heart Study and MESA CV risk indices (Figure 3, top panels).
None of the traditional time (AVNN, SDNNIDX, rMSSD, pNN50) and frequency domain (HF, LF/HF) HRV variables was significantly associated with incident events. DFA α 1 was also not associated with incident events.  Model 1: unadjusted. Model 2: adjusted for the traditional risk factors: age, sex, systolic blood pressure, total cholesterol, HDL cholesterol, current smoking status, hypertension medication (including beta blockers, calcium-channel blockers, angiotensin antagonists, and angiotensin antagonists plus diuretics), diabetes and lipid lowering medication. Values presented are standardized hazard ratios (HRs), 95% confidence intervals (95% CI), Harrell's C statistic (C-index) and the p-value for the likelihood ratio test (LR-test) of the null hypothesis that the addition of a dynamical measure (fragmentation, traditional HRV or DFA α 1 metric) to a model with the traditional risk factors did not improve the fit of the data. The numbers of participants/events in the analyses of the models 1 and 2 were 1771/72 and 1702/71, respectively. Abbreviations: Q 1 , first quartile; Q 3 , third quartile; HRV, heart rate variability; CVE, cardiovascular event; HDL, high density lipoprotein; PIP, percentage of inflection points; W 0 , W 1 , W 2 , W 3 , percentage of words with 0, 1, 2, and 3 inflection points, respectively; AVNN, average value of the NN intervals; SDNNIDX, mean of the standard deviations of NN intervals in all 5-min segments; rMSSD, root mean square of the successive differences; pNN50, percentage of differences between successive NN intervals above 50 ms; HF, high frequency spectral power; LF/HF, ratio of low to high frequency power; DFA α 1 , detrended fluctuation analysis short-term fractal exponent. Boldface is used to highlight p-values < 0.05.

Analyses of Risk of Incident CVEs Adjusted for Traditional Risk Factors
The models with each of the fragmentation and traditional HRV variables were adjusted for standard CV risk factors: age, sex, systolic blood pressure, total cholesterol, HDL cholesterol, current smoking status, hypertension medication, diabetes and lipid lowering medication. In these analyses ( Table 2: Model 2), all fragmentation indices remained significantly associated with incident CVEs. For example, a one-SD increase in PIP was associated with a 43% (15-78%) increase in the rate of incident CVEs. In addition, all fragmentation indices, with the exception of W 3 calculated from the RR interval time series, added significant value to the model with only the risk factors.
In contrast, none of the traditional time (AVNN, SDNNIDX, rMSSD, pNN50) and frequency domain (HF, LF/HF) HRV variables was significantly associated with incident events. DFA α 1 was also not associated with incident events. Age and sex were the only covariates significantly associated with incident CVEs in the adjusted models that included a fragmentation index. Specifically, for a model with PIP, the standardized hazard ratios for each of the covariates were: i) age (SD = 9.1 years), 1.33 (95% CI: 1.03-1.73), p = 0.029; ii) sex, 1. The results did not qualitatively change after further adjusting the analyses for each of the following variables, singly or in combination: race, body mass index, waist circumference, diastolic blood pressure, seated heart rate, use of hypoglycemic agents, total sleep time, sleep efficiency, the apnea-hypopnea index, measures of coronary artery calcification (McClelland et al., 2015) (total Agatston calcium score) and carotid plaque (Gepner et al., 2015(Gepner et al., , 2017.
For all sleep stages, including awake during sleep, and awake before sleep onset and/or sleep termination, the degree of HRF, measured by PIP, was positively associated with risk of incident

Analyses of Risk of Incident CVEs Adjusted for the Framingham and MESA CV Risk Indices
In general, the fragmentation indices remained significantly associated with the risk of CVEs in models adjusted for the Framingham and the MESA CV risk indices (Table 3). Specifically, increased fragmentation, that is, higher PIP, lower percentages of fluent words W 0 and W 1 , and higher percentages of fragmented words W 2 and W 3 , were significantly associated with increased risk of events. Neither the traditional HRV measures nor DFA α 1 showed significant associations with incident CVEs.
The risk indices in each of these models were also significantly associated with incident CVEs. Specifically, a one-SD increase in the Framingham and in the MESA risk indices, was associated with 80% (95% CI: 43-125%), and 55% (33-81%) increase in the rate of adverse CVEs, respectively. Harrell's C statistic was 0.666 and 0.678 for the Framingham and MESA risk indices, respectively. Overall the best model, with a Harrell's C statistic Model 3: adjusted for the Framingham CV risk index (D'Agostino et al., 2008). Model 4: adjusted for the MESA CV risk index. Values presented are standardized hazard ratios (HRs), 95% confidence intervals (95% CI), Harrell's C statistics (C-index) and the p-value for the likelihood ratio test (LR-test) of the null hypothesis that the addition of an HR dynamical metric (fragmentation, traditional HRV or DFA α 1 ) to a model with the Framingham or MESA risk indices did not improve the fit of the data. The numbers of participants/events in the analyses of models 3 and 4 were 1767/72 and 1672/71, respectively. Abbreviations: HRV, heart rate variability; CVE, cardiovascular event; CV, cardiovascular; PIP, percentage of inflection points; W 0 , W 1 , W 2 , W 3 , percentage of words with 0, 1, 2, and 3 inflection points; AVNN, average value of the normal-to-normal sinus (NN) intervals; SDNNIDX, mean of the standard deviations of NN intervals in all 5-min segments; rMSSD, root mean square of the successive differences; pNN50, percentage of differences between successive NN intervals above 50 ms; HF, high frequency spectral power; LF/HF, ratio of low to high frequency power; DFA α 1 , detrended fluctuation analysis short-term fractal exponent. Boldface is used to highlight p-values < 0.05. of 0.703, was the one that combined the word group W 1 derived from RR intervals, with the MESA risk index. Models incorporating the Framingham risk index and any of the fragmentation measures, except W 2 derived from the NN interval time series, performed significantly better than the Framingham index itself. Similarly, all models that included any of the fragmentation measures in addition to the MESA risk index performed significantly better than the MESA index itself.
The traditional HRV variables and α 1 were not significantly associated with risk of incident CVEs either in unadjusted or adjusted models. Adding one of these indices to a model with a fragmentation index did not improve its performance.

Analyses of Risk of CV Death: Unadjusted and Adjusted for the Framingham and MESA Risk Indices
Higher PIP and lower percentage of W 1 words were significantly associated with increased risk of CV death in unadjusted analyses as well as in analyses adjusted for the Framingham and the MESA CV risk indices (Table 4, Figure 3). Specifically, a one-SD (∼ 7%) increase in PIP derived from the analysis of NN interval time series was associated with an increase in the rate of CV death of 89% (95% CI: 34-168%) in unadjusted models and of 65% (15-136%) and 67% (19-136%) in models adjusted for the Framingham and the MESA risk indices, respectively. A one-SD (∼ 11%) increase in the percentage of W 1 ("fluent" or least fragmented) words, also derived from the analysis of NN interval time series, was associated with a decrease in the rate of CV death of 59% (37-75%) in unadjusted models and of 52% (21-71%) and 55% (25-72%) in models adjusted for Framingham and the MESA risk indices, respectively. Similar results were obtained from the analyses of the RR interval time series. Further adjusting the models by prevalent CVD did not change the significance of the associations between the fragmentation metrics and risk of CV death.
Lower percentages of fluent words W 0 were associated with increased risk of CV death in all models. However, these associations were weaker than those with word W 1 . The percentages of word W 2 and W 3 , the most fragmented, were positively associated with the risk of CV death in all models. However, the significance of the associations depended on the particular model (Table 4). A one-SD increase in the Framingham (∼9%) and in the MESA risk (∼6%) indices was associated with 137% (95% CI: 48-278%) and 62% (35-96%) increase in the rate of CV death, respectively. Harrell's C statistic for the former and latter models was 0.749 and 0.797, respectively. Both variables, PIP and W 1 , calculated from NN and RR interval time series, added significant information to models with either of these risk indices. The results for word groups W 0 , W 2 , and W 3 depended on the particular model. Overall,  (HRs), 95% confidence intervals (95% CI), Harrell's C statistics (C-index) and the p-value for the likelihood ratio test (LR-test) of the null hypothesis that the addition of a dynamical measure (fragmentation, traditional HRV or DFA α 1 metric) to a model with one of the risk indices did not improve the fit of the data. Abbreviations: HRV, heart rate variability; CV, cardiovascular; PIP, percentage of inflection points; W 0 , W 1 , W 2 , W 3 , percentage of words with 0, 1, 2, and 3 inflection points; AVNN, average value of the normal-to-normal sinus (NN) intervals; SDNNIDX, mean of the standard deviations of NN intervals in all 5-min segments; rMSSD, root mean square of the successive differences; pNN50, percentage of differences between successive NN intervals above 50 ms; HF, high frequency spectral power; LF/HF, ratio of low to high frequency power; DFA α 1 , detrended fluctuation analysis short-term fractal exponent. Boldface is used to highlight p-values < 0.05.
the best model, with a Harrell's C statistic of 0.838, was the one that combined the word W 1 , derived from RR intervals, with the MESA risk index. The traditional HRV variables and α 1 were not significantly associated with risk of CV death either in unadjusted or adjusted models. Adding one of these indices to a model with a fragmentation index did not improve its performance.

Relationship of HRF With Short-Term Traditional HRV Indices
Nonlinear (U-shape) relationships were found between fragmentation indices and traditional HRV measures of shortterm variability. Figure 4 shows one representative example, the relationship between PIP and ln(rMSSD). For the first three quartiles of rMSSD values, (i.e., for rMSSD values below the 75th percentile of rMSSD, specifically, ln(rMSSD) < 3.7 ms), the degree of fragmentation and the amount of short-term variability were inversely correlated. In the upper quartile of rMSSD values, however, the degree of fragmentation and the amount of short-term variability were positively associated. Qualitatively similar results were found for pNN50 and HF power.

DISCUSSION
The present investigation was designed to test the association of quantitative measures of HRF, a newly defined property of shortterm sino-atrial rhythm dynamics, with adverse CV outcomes in MESA, a large ongoing multicenter study of individuals recruited from the general community. The key findings of this study were that: (1) increased HRF was significantly associated with risk of incident CVEs and CV mortality; (2) measures of fragmentation added value to Framingham and MESA risk prediction indices; and (3) traditional metrics of short-term HRV as well as a nonlinear index (DFA α 1 ) were not associated with incident CVEs or CV death. The development of the concept of fragmentation and its quantitative metrics was motivated in part by an apparent paradox in the results of traditional time and frequency domain analyses of a number of studies (de Bruyne et al., 1999;Stein et al., 2005;Huikuri and Stein, 2012;Drawz et al., 2013;Costa et al., 2017a;Raman et al., 2017). While the mechanisms of the relatively slow (i.e., below the respiratory frequency) variations in HR are attributable to complex interactions between the parasympathetic and the sympathetic branches of the autonomic nervous system, faster variations, in the range of 0.15 − 0.40 Hz, are mainly attributed to vagal tone modulation (HRV, 1996). Therefore, short-term (high frequency) measures of HR dynamics, such as rMSSD, pNN50, and HF power, are typically interpreted as surrogate measures of cardiac vagal tone. In contexts where cardiac vagal tone modulation is known to be diminished, for example, with advanced aging and established CVD, these "vagal" measures are expected to be lower. In fact, a monotonic decrease in high frequency variability with increasing age has been reported in multiple cross-sectional studies of ostensibly healthy adults (Pikkujamsa et al., 1999;Bonnemeier et al., 2003;Costa et al., 2017a). Furthermore, the association between extremely low variability and adverse outcomes is welldocumented (HRV, 1996;Huikuri and Stein, 2012;O'Neal et al., 2016).
However, in certain cases (de Bruyne et al., 1999;Stein et al., 2008;Almeida-Santos et al., 2016;Raman et al., 2017;Wdowczyk et al., 2018) a paradoxical increase in short-term HR fluctuations was observed in contexts where reduced vagal tone would have been expected based on age and/or advanced heart disease. In the present study, we observed a U-shaped relationship between traditional short-term HRV measures and cross-sectional age (Figure 2). From approximately ages 45-65 years the amount of short-term variability decreased. Subsequently, variability increased despite the well-known decrease in cardiac vagal tone modulation with advancing age. These results provide further evidence that in cohorts of middle-aged to elderly individuals, such as MESA, traditional HRV indices may fail to reflect accurately changes in cardiac vagal tone.
We introduced fragmentation analysis to quantify short-term HR dynamics in such types of cohorts. The term "fragmented heart rate" refers to rhythms in which HR acceleration sign changes at a frequency higher than that attributable to healthy vagal tone modulation of the SAN. These anomalous rhythms include but are not limited to classic sinus alternans and its variants. Due to their anti-correlated structure, i.e., the fact that an increase in HR is likely followed by a decrease and viceversa, fragmented sinus rhythms are not random/erratic. In fact, they exhibit a much higher degree of predictability/regularity (Figure 1, bottom panels) than random rhythms. The term "fragmented" was also chosen to help convey the putative pathophysiologic concept of regulatory network disintegration and/or breakdown of neuroautonomic coupling between heart rate and respiration.
An intuitive measure of fragmentation is the percentage of changes in HR acceleration sign, that is, PIP, in NN (or RR) time series. Most recently, a symbolic dynamical approach was introduced (Costa et al., 2017b) that quantifies the frequency of occurrence of different patterns of fluctuations, from least fragmented (most fluent) to most fragmented. These fragmentation indices were originally tested in studies of publicly available databases from the Rochester THEW archives. Of note, if the amplitude of the fluctuations is low (e.g., 80 ms), fragmentation is unlikely to be detected in clinical readings of short (typically 10 s) and long (Holter) ECG recordings. However, as shown in Figure 1, visual detection of HRF is facilitated by inspection of HR time series.
The origins of HRF remain speculative. Possible pathophysiologic mechanisms include increased automaticity in or proximal to the SAN, exit block in the SAN area, modulated sinus/atrial parasystole, abrupt pacemaker shifts in the SAN (Boyett et al., 2000;Kodama et al., 2004), beat-to-beat changes due to perturbations in atrial stretch receptors (Costa et al., 2017a,b) or alterations in membrane and cellular pacemaker clocks (Lakatta et al., 2010). These electrophysiologic perturbations, in turn, may be related to underlying atrial (Sosnowski and Petelenz, 1995;Roberts-Thomson et al., 2007;Goette et al., 2017) or ventricular disease. Systemic and local factors that may also contribute to pathophysiologic dysregulation of SAN dynamics include inflammation, degeneration, fibrosis and calcification (Costa et al., 2017a). Future experimental and mathematical modeling studies will hopefully shed light on the putative links between these and other mechanisms and fragmentation. Possible genomic associations with HRF remain to be explored.
Based on the analyses of the THEW databases (Costa et al., 2017a,b), we hypothesized that increased HRF might be a biomarker of increased risk of incident CVEs and CV death. To explore these hypotheses, we analyzed HR dynamics from a subset of the participants in the MESA. This national study is one of the largest prospective investigations designed to track meticulously the course of CVD in an ethnically diverse population free of overt clinical CVD at study entry. Two types of ECG recordings with detailed follow-up data were available at the time of this analysis: traditional 10-s ECGs and the ECG channel of the PSG studies. We chose to examine the latter given the non-stationary nature of HR dynamics, [which implies that statistical time series analysis tools are most reliable when applied to "long" recordings (HRV, 1996)], and the fact that previous studies (Costa et al., 2017a,b) have shown that the discriminatory power of HRF was generally comparable during awake and sleep periods. In MESA, we confirmed and extended this initial observation. Specifically, we found that the positive association between HRF and incident CVEs (in fully adjusted analyses) could be detected for each of the sleep stages and for awake periods before sleep onset and/or after sleep termination.
A concomitant finding was the absence of associations between the most commonly used traditional short-term HRV measures and incident CVEs and CV death, in unadjusted and adjusted models. These results are not as surprising as they might appear at first glance. First, the U-shaped relationship between traditional short-term HRV measures and the participants' age (Figure 2) was indicative that such measures would be of limited utility in this cohort. Second, as previously mentioned, traditional measures of HRV, in contrast to fragmentation measures, also failed to discriminate patients with CAD from ostensibly healthy subjects in databases provided by the University of Rochester (Costa et al., 2017a,b). Third, HRF, by increasing variability not ascribable to physiologic vagal tone modulation may confound the results of traditional HRV. The nonlinear (U-shaped) relationship between HRF and short-term variability (Figure 4) supports this conjecture. In fact, the subgroups of participants with the lowest and highest amounts of variability, which, in the conventional HRV framework, would be presumed to have the highest and lowest risk of adverse events, respectively, both showed increased HRF. These findings are consistent with the reports of Stein et al. (Stein, 2002;Stein et al., 2008) and others (de Bruyne et al., 1999;Huikuri and Stein, 2012;Drawz et al., 2013).
Of note, fragmentation and traditional HRV indices differ in the following major way. By construction, fragmentation indices do not mathematically depend on mean HR and/or the amplitude of its fluctuations. These salient attributes derive from the fact that accelerations/decelerations are defined as increments/decrements in HR of any magnitude. In contrast, by definition, short-term HRV indices quantify information that is encoded in the amplitude of the fluctuations. As previously mentioned rMSSD, pNN50, and HF power were not associated with risk of incident CVEs and CV mortality. The other widely used HRV metrics, AVNN, SDNNIDX, and LF/HF were also not associated with risk of incident CVEs and CV mortality. Furthermore, none of the traditional indices improved the performance of models that included a fragmentation index.
In this study the short-term DFA index, α 1 , decreased with the participants' age. In addition, α 1 was lower in those with adverse events (incident CVEs and CV death) than in those without. However, the associations were not statistically significant in any of the models.
Of potential basic and translational importance is the fact that we used both the NN and RR time series in fragmentation analyses. The NN series were employed using expert edited time series from THEW and MESA to insure that the fragmentation was likely related to beats originating in or near to SAN, therefore not distinguishable from sinus beats, at least from the single lead provided. The RR time series were used to demonstrate that fragmentation analysis, not relying on detailed beat annotation, had comparable (or even superior) discriminatory power to that employing NN time series, substantially facilitating the development of automatable analyses.
Finally, it is worth emphasizing that the Framingham and the MESA indices are composite measures incorporating information related to demographics (age, sex, race), lifestyle (smoking status), vital signs (blood pressure) and blood analytes (lipids, glucose). In contrast, each HRF index is a dynamical measure reflecting the frequency of the changes in HR acceleration sign. How can such a single metric based on a continuous ECG keep "pace" with these other multivariable risk stratification tools? The answer may relate in part to the fact that HRF indices are dynamical measures, not static probes. In contrast, blood pressure, cholesterol, glucose and other common biomarkers are single time point readouts ("snapshot"). Thus, they provide limited information on the dynamics of the underlying control mechanisms.
More generally, HRF metrics belong to the new class of dynamical probes (Cysarz et al., 2000;Goldberger, 2006;Costa et al., 2014;Porta et al., 2015;Makowiec et al., 2015;Hoyer et al., 2017;Zhang et al., 2017) that can report, in real-time, on salient aspects of integrative, multiscale, regulatory systems and of their breakdown with aging and disease. The use of these probes may enhance the clinical utility of traditional risk assessment tools (Tables 3, 4) and of other emerging technologies, such as genomic profiling. In furtherance of the goals of precision medicine, the dynamical property of HRF may also constitute a novel "target" for therapeutic interventions.

CONCLUSION
HRF, a newly defined manifestation of anomalous shortterm sino-atrial variability, is associated with increased risk of cardiac adverse events and cardiac mortality in MESA. The measures, readily computable from long, continuous ECGs, added value to the canonical Framingham and MESA CV risk indices. Furthermore, fragmentation measures outperformed conventional HRV measures and the short-term DFA α 1 index. Future studies are needed to confirm the utility of HRF measures in risk stratification, and in prediction of cardiomyopathies and atrial fibrillation.

AUTHOR CONTRIBUTIONS
The idea for this study was developed by MC and AG. Data acquisition was the responsibility of SR. The analysis of the data was performed by MC. RD and SH directed the statistical analysis. All co-authors participated in the interpretation of the findings. MC and AG wrote the manuscript in collaboration with the other co-authors.